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The  role  of  the  Maneuvering  Target  Statistical  Tracker  (MTST),  a  Kalman  filter  tracking 
algorithm  based  on  the  Integrated  Ornstein-Uhlenbeck  (IOU)  motion  process,  in  the 
Antisubmarine  Warfare  System  Evaluation  Tool  (ASSET)  is  examined  and  its  operation 
described.  ASSET  is  a  campaign  simulation  which  models  open-ocean  ASW  scenarios  featuring 
prosecution  of  hostile  submarines  by  friendly  submarines  and  aircraft  based  on  cues  provided 
by  data  fusion  centers.  The  heart  of  each  data  fusion  center  is  an  MTST  which  integrates  new 
contact  information  into  tracks.  Comparing  the  level  of  sophistication  of  the  tracking  algorithm 
to  that  of  the  contact  data  provided  to  it,  a  number  of  simplifications  are  proposed.  These 
include  using  reduced  complexity  IOU  prediction  and  Kalman  filter  equations;  the  use  of  pre- 
processed  variance  data  together  with  the  true  position  of  targets  to  estimate,  rather  than 
explicitly  calculate,  updated  track  states;  and  limiting  contact  processing  based  on  information 
content.  Results  indicate  a  good  simulation  of  tracker  output  is  produced  using  a  greatly 
simplified  algorithm.  This  technique  can  be  generalized  to  other  types  of  simulations  involving 
target  tracking. 
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I.     AN  INTRODUCTION  TO  ASSET 

The  Antisubmarine  Warfare  Systems  Evaluation  Tool  (ASSET)  is  a  campaign-level 
simulation  which  models  open  ocean  scenarios  involving  submarines,  maritime  patrol 
aircraft  (MP A),  shore-based  command  and  data-fusion  centers  and  a  wide  variety  of 
passive  acoustic  sensors.  Active  acoustic  and  non-acoustic  sensors  are  only  modeled 
using  a  simple  area  sensor  model  which  has  a  fixed  probability  of  detection  and  false 
alarm  rate.  The  simulation  scenario  is  specified  by  a  user-supplied  "architecture"  which 
determines  all  facets  of  environment,  command  control,  sensor  interaction  and  platform 
maneuver.  This  structure  is  input  to  a  desktop  computer/workstation  (the  Apple 
Macintosh  II  in  this  version)  through  a  series  of  user-friendly  windows,  each  dealing  with 
a  specific  topic.  A  particular  configuration  of  sensors,  data  fusion  centers, 
communication  nodes  and  tactical  platforms  interact  in  a  particular  geographic  location 
against  an  analogous  enemy  force  structure. 

The  scenario  is  repeated  as  a  Monte  Carlo  simulation  to  produce  statistically 
meaningful  measures  of  effectiveness  (MOE's).  Output  data  regarding  the  detection, 
localization  and  prosecution  of  enemy  submarines  can  then  provide  a  quantitative  basis 
for  decisions  regarding  ASW  Master  Plans,  Top  Level  Warfare  Requirements,  and  a 
variety  of  emerging  technology  assessments  and  system  appraisals  [Ref.  l:p.  iv].  These 
results  can  also  be  useful  in  conducting  quantitative  experimentation  with  regional  force 


levels,  force  compositions  and  commitment  strategies;  command,  control, 
communications  and  intelligence  (C3I)  networks;  implications  of  foreign  technology 
advances;  and  fleet  exercise  planning.  The  current  version  of  ASSET  (version  1 .0)  limits 
the  scope  of  these  scenarios  to  open  ocean  search  and  prosecution  of  hostile  submarines 
by  cueing  friendly  submarines  and  MPA  with  wide-area  sensors,  and  modeling  the 
supporting  C3I  networks.  The  modular  nature  of  the  object-oriented  structure  of  ASSET 
makes  expanding  the  scope  of  the  simulation  possible  [Ref.  2:p.  1-2].  As  the 
construction  and  operating  costs  of  naval  platforms  spiral  upwards,  a  simulation  tool  of 
this  kind  will  be  increasingly  important  to  intelligently  manage  resources  vis-a-vis  a 
rapidly  evolving  spectrum  of  possible  threats. 

A.      OBJECTIVES  AND  RATIONALE 

The  purpose  of  this  thesis  is  three-fold: 

•  to  provide  a  description  of  the  target  tracking  algorithm  used  by  ASSET  and 
how  it  relates  to: 

1.  real- world  tracking. 

2.  the  ASSET  simulation  as  a  whole. 

•  to  examine  the  mathematical  development  of  the  tracking  algorithm. 

•  to  examine  possible  modifications  to  the  existing  Kalman  filter  tracking 
algorithm  to  better  match  the  data  input  to  it  and  increase  computational 
efficiency. 


The  use  of  full  Kalman  filter  equations  for  elliptical  areas  of  uncertainty 
(AOU's)  is  not  necessary  given  the  circular  AOU's  that  ASSET'S  sensor  model 
produces.   Making  use  of  the  fact  that  the  computer  knows  the  ground  truth  locations 
of  the  platforms,  and  the  limited  data  input  to  the  filter,  it  will  be  shown  that  it  is 
possible  to  approximate  appropriately  distributed  filtered  state  positions  using 
simplified  forms  of  the  Kalman  filter  matrix  computations.   The  situations  when  the 
tracker  algorithm  is  used  can  also  be  limited  based  on  the  information  content  of  the 
contact.   Together,  these  modifications  yield  a  significant  reduction  in  the 
computational  complexity  of  the  tracking  algorithm,  which  is  the  principle  goal  of  this 
thesis. 

The  automatic  correlator  tracker  (ACT)  implemented  in  ASSET  is  a  stand-alone 
module  based  on  the  Ocean  Surveillance  Information  System  (OSIS)  baseline  upgrade 
single  hypothesis,  multiple  target,  Kalman  filter-based,  correlation  and  tracking 
algorithm  [Ref.  3:p.  4].  This  set  of  routines  contains  a  complete  representation  of  the 
OSIS  baseline  upgrade  automatic  correlator  tracker's  (OBU-ACT)  functions,  with  the 
exception  of  LINK,  the  utility  that  evaluates  track  sets  which  are  potentially  legs  of  a 
single  target's  track;  and  EQUATE,  which  makes  track  associations  for  contacts  pre- 
correlated  to  a  particular  platform  based  on  acoustic  or  electromagnetic  signature. 

The  OBU-ACT  package  is  designed  to  process  contact  information  derived  from 
position-only,  bearing-only,  and  position  and  velocity  reports.   The  covariance 
matrices  associated  with  these  reports  can  be  interpreted  as  an  elliptical  representation 
of  the  error  in  position  and  velocity.  The  capability  to  process  line  of  bearing 


contacts  is  not  available  in  ASSET  (1.0).  The  method  used  to  model  sensor  systems 
in  ASSET  1.0  approximates  the  elliptical  errors  with  circular  ones.   This  reduces  the 
complexity  of  the  covariance  matrices  considerably.   Taking  advantage  of  the  large 
number  of  zero  elements  and  repeated  values  in  the  matrices  actually  constructed  by 
ASSET,  the  matrix  equations  involved  in  the  tracking  algorithm  can  be  reduced  to 
equivalent  scalar  ones  by  eliminating  all  the  zero-multiplied  terms. 

Central  to  this  algorithm  is  the  Integrated  Ornstein-Uhlenbeck  (IOU)  process 
which  models  target  motion  in  the  Maneuvering  Target  Statistical  Tracker  (MTST). 
The  MTST  uses  this  model,  together  with  contact  data  corrupted  by  noise,  to  estimate 
the  size  and  position  of  the  Submarine  Probability  Area  (SPA)  which  has  an  86 
percent  probability  of  containing  the  target's  true  position.   The  major  modification  to 
the  algorithm  proposed  is  based  on  assuming  that  the  IOU  prediction  position 
distribution,  in  like  fashion  to  the  contact  position  distribution,  is  centered  on  the 
target's  true  position.  This  results  in  a  distribution  of  SPA  centers  that  is  centered  on 
the  targets  true  position  also.   Thus,  only  the  SPA  variance  is  calculated,  not  its 
center  position,  which  is  drawn  randomly  from  the  estimated  distribution  about  the 
true  target  position.   Since  contact  data  does  not  play  a  role  in  computing  these 
variances,  they  may  be  calculated  prior  to  the  start  of  the  simulation.   A  table  of 
variances  for  the  sensor  AOU  values  defined  for  the  scenario  and  an  appropriate 
range  of  contact  report  interarrival  times  can  be  constructed  and  referred  to  as  needed 
during  simulation  execution  to  minimize  the  time  spent  processing  contact  data. 


B.      FUNDAMENTALS  OF  OPERATION. 

The  basic  structure  of  ASSET  was  adapted  from  COAST,  the  Common  LISP 
Architectural  Study  Tool.   Common  LISP  is  an  object-oriented  programming  language 
used  for  rapid  prototyping  and  artificial  intelligence  work.   A  general  overview  of  this 
structure  is  presented  to  facilitate  understanding  of  how  the  COAST/ASSET  system  is 
organized  and  the  critical  role  of  the  tracker  to  simulation  operation. 

1.     Objects  and  Object  Interactions. 

Object-oriented  programs  combine  data  with  associated  procedures  to  form  a 
hierarchical  network  of  self-contained  modules  known  as  objects.   These  objects  can 
then  be  invoked  by  each  other  according  to  the  program  methodology.   Program 
objects  are  independent  and  can  be  modified  without  affecting  other  objects  they 
communicate  with.   The  data  that  is  communicated  will  change,  but  the  relationship 
between  the  objects  does  not  necessarily  have  to.  Data  can  be  evaluated  internally  by 
an  object  without  affecting  the  analogous  data  in  a  separate  object. 

The  principle  of  instancing  allows  a  single  object  definition  to  create  any 
number  of  structures  within  the  program,  all  obeying  the  same  definition.  This 
allows  a  single  definition  of  a  submarine  object,  for  example,  to  be  given  many 
different  sets  of  parameters  each  representing  a  different  class  of  submarine.   Other 
objects  may  be  designated  to  impart  their  characteristics  to  the  new  submarine  object. 
This  process  of  inheritance  allows  the  properties  of  another  object,  an  acoustic  sensor 
object  for  example,  to  be  created  separately  from  the  submarine  object.  The 


submarine  object  can  then  be  designated  to  be  a  kind  of  acoustic  sensor  and  it  will 
inherit  the  properties  of  the  acoustic  sensor  object  in  addition  to  its  own.   Another 
object  such  as  a  SURTASS  towed-array  surveillance  ship  object  can  also  inherit  the 
properties  of  the  acoustic  sensor  object,  with  different  parameters,  without  affecting 
the  parameters  of  the  acoustic  sensor  inherited  by  the  submarine  object. 

Complex  classes  of  platforms  can  be  formed  by  multiple  inheritance  of  the 
properties  of  basic  building  block  objects.  These  are  then  used  to  simulate  any 
number  of  platforms  of  that  class,  each  operating  independently  within  the  simulation. 
The  computer  memory  available  to  store  each  of  these  instances  of  the  object  becomes 
the  limiting  factor  to  the  complexity  of  the  scenario  to  be  examined.  The  organization 
of  objects  in  ASSET  consist  of  several  major  groups: 


•  Graphical  environment  objects  which  implement  the  standard  Macintosh 
graphical  user  interface  to  provide  input/output  by  means  of  a  mouse  based 
"point  and  click"  metaphor  featuring  pull-down  menus, dialog  boxes  with  labeled 
data  regions  and  radio  buttons,  and  windowed  presentation  of  multiple 
information  sources  simultaneously. 

•  Region  management  objects  which  construct  and  manage  distinct  regions  of 
varying  environmental  properties,  command  responsibility  or  platform  patrol 
assignment. 

•  Event  management  objects  which  queue  all  simulation  events  in  time  order 
sequence  and  parcel  them  out  to  the  appropriate  resolution  objects. 

•  Command  objects  which  perform  resource  allocation,  assigning  available  assets 
to  contact  cues  based  on  either  time  to  station  (MP A)  or  area  of  uncertainty  size 
(submarines). 

•  Automatic  Correlator  Tracker  (ACT)  objects  which  manage  the  grouping  and 
processing  of  contact  reports  into  target  tracks  creating  a  tactical  picture 
consisting  of  combinations  of  true  and  false  contacts. 


•  Tactical  platform  objects  which  simulate  the  behavior  of  their  real  world 
counterparts.   These  are  dynamic  associations  of  component  objects  which  may 
change  in  response  to  simulation  events.   A  submarine  inheriting  the  properties 
of  a  Patroller  object  may  switch  to  an  Intercepter  after  making  a  detection. 


The  program  flow  from  object  to  object  is  not  sequential,  but  event  driven. 
Once  started,  the  simulation  clock  advances  time  relative  to  the  simulation.   Sensor 
objects  begin  glimpsing  and  motion  platforms  begin  moving.   Possible  detection 
events  are  sent  to  the  event  manager  for  proper  sequencing  and  resolution.   Detections 
are  communicated  up  the  designated  chain  of  command  to  data-fusion  centers  where 
an  instance  of  the  ACT  processes  them  and  integrates  them  into  the  tactical  picture. 
This  may  cause  a  command  object  to  direct  (or  redirect)  assets  in  response.   This 
process  continues  until  the  allotted  time  for  the  scenario  expires.   The  MOE  statistics 
for  that  run  are  added  to  the  MOE  file  and  the  next  iteration  of  the  scenario 
commences.   When  the  desired  number  of  iterations  are  complete  the  compiled 
summary  of  MOE  statistics  is  analyzed  to  interpret  the  results  of  the  simulation. 

2.     ASW  System  Architecture  Evaluation. 

In  order  to  evaluate  a  desired  architecture  it  must  be  conceptualized  in  great 
detail.  The  classic  computer  maxim  "garbage  in,  garbage  out"  is  especially  true  of 
ASSET  where  a  single  bad  parameter  value  can  render  the  results  meaningless. 
While  the  user-friendly  interface  facilitates  the  mechanical  process  of  inputing  data, 
the  abstract  nature  of  the  required  parameters  make  it  necessary  to  assume  values  of 
dubious  validity  at  times.    General  areas  which  require  quantified  data  include: 


Communication  connectivities  which  represent  how  detection  reports  are  passed 
from  object  to  object  and  what  delays  are  involved  at  each  node. 

Command  organization  including  geographic  regions  of  responsibility. 

Environmental  data  for  acoustic  propagation  loss  as  a  function  of  range  and 
frequency  as  well  as  ambient  acoustic  noise  in  each  distinct  environmental 
region. 

Motion  plans  representing  the  operating  areas  and  interconnecting  tracks  which 
will  govern  tactical  platform  movement. 

Umpire  parameters  such  as  the  kill  probabilities  for  a  given  submarine  class 
against  each  possible  target  class  both  for  the  case  where  it  detects  the  enemy 
first  and  the  case  where  it  is  detected  by  the  enemy  first. 


The  individual  platforms  require  complex  definition  as  well.   A  submarine 
object,  for  example,  requires  parameters  for: 

•  breakoff  speed  between  "fast"  and  "slow"  acoustic  behavior. 

•  self-noise  at  slow  and  fast  speed. 

•  directivity  index  and  recognition  differential  (for  Sonar  Equation  computations). 

•  signature  frequency  emitted  and  intensity  at  slow  and  fast  speed. 

•  movement  speeds  when  patrolling  and  intercepting  and  detailed  motion  plan  to 
be  followed. 

•  weapons  loadout  and  level  at  which  the  sub  will  abort  its  mission  to  rearm. 

•  whether  the  sub  will  transmit  the  detections  it  makes  and  risk  detection  itself  or 
not  report  any  detections. 

•  the  interval  at  which  it  copies  the  submarine  broadcast  for  orders  to  investigate 
cues.  [Ref.  2:p.  2-37] 
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The  submarine's  sensor  requires  a  detailed  set  of  parameters  of  its  own  which 
will  be  used  to  construct  the  contact  report  that  is  communicated  to  the  ACT 
representing  the  data  fusion  center  associated  with  the  submarine's  command  object: 


•  position  error  which  represents  the  2  standard  deviation  radius  of  a  circular 
normal  distribution  of  contact  errors  about  the  actual  target  position. 

•  course  and  speed  errors  which  are  uniform  about  the  true  value. 

•  Target  Motion  Analysis  (TMA)  delay  representing  the  time  required  to  acquire 
course  and  speed  information. 

•  false  alarm  rate. 

•  determination  whether,  on  the  average  the  sensor  is  reliable  enough  for  the 
fusion  center  to  automatically  start  a  new  track  based  on  a  single  contact  report. 


The  complexity  of  the  data  involved  in  constructing  an  architecture  for 
evaluation  varies  considerably,  as  can  be  seen  from  these  examples.   Once  all  this 
information  is  input,  it  can  be  easily  double  checked  by  accessing  the  appropriate 
windows  again.   When  the  user  is  sure  all  data  is  input  correctly,  the  simulation  can 
be  run  with  the  platform  graphics  turned  on  to  ensure  it  progresses  properly.   These 
graphics  can  then  be  turned  off  to  run  multiple  iterations  faster. 

When  the  desired  number  of  iterations  is  completed  a  dialog  box  opens  asking 
if  the  MOE's  are  to  be  saved.   Once  saved  to  disk,  they  can  be  opened  and  examined. 
They  include: 


•  Attrition  of  submarines  and  MPA  assets  for  each  side  as  a  fraction  of  original 
force. 

•  Number  of  times  submarines  approached  to  within  a  critical  range  (assigned  by 
the  user)  of  enemy  surface  formations  (which  serve  only  as  targets  in  ASSET 
1.0,  they  have  no  combat  or  detection  capability). 

•  Tracker  statistics  involving  the  fusion  delay  involved  in  constructing  target 
tracks. 


These  statistics  along  with  weapon  expenditures  allow  the  completed  simulation 
to  be  compared  to  other  runs  to  observe  the  effect  of  a  particular  parameter  changing, 
or  how  close  a  particular  MOE  comes  to  a  goal  value.   The  methodology  of  ASSET 
must  be  taken  into  account  as  the  computer  commanders  follow  simple  resource 
allocation  rules  which  do  not  necessarily  reflect  the  tactical  priorities  the  user  desires. 
A  higher  priority  for  prosecution  may  well  be  given  to  a  submarine  returning  to  base 
than  one  ten  miles  from  a  surface  group,  due  to  the  resource  allocation  process' 
emphasis  on  detection  rather  than  on  protection.   While  creativity  is  needed  in  the 
design,  implementation  and  interpretation  of  an  ASSET  scenario,  significant  insights 
into  the  conduct  of  ASW  campaigns  can  be  gleaned  from  it. 

3.     Correlation,  Tracking  and  Data  Fusion. 

ASSET  is  designed  to  simulate  the  data  fusion  centers  where  contact  reports 
from  a  wide  variety  of  sensors  and  platforms  are  centrally  processed  to  create  a 
tactical  picture  for  the  region  of  interest.   These  data  fusion  centers  are  simulated  by 
an  instance  of  the  OBTJ-ACT  module  with  appropriate  communication  connectivities. 
The  result  of  correlating  and  processing  the  streams  of  contact  information  arriving  at 
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a  fusion  center  forms  the  basis  for  the  allocation  and  cuing  of  assets.   Thus  it  is 
critical  that  the  ACT  is  not  spinning  its  wheels  performing  unnecessary  computations. 

The  implementation  of  OBU-ACT  used  in  ASSET  is  organized  into  four  main 
functional  areas: 


•  generating  a  new  track  from  one  or  more  contact  reports  (the  START  and 
CLUSTER  modules). 

•  associating  a  contact  report  to  an  existing  track  (the  INPUT,  START,  and 
ASSOC  modules). 

•  simulating  the  intervention  of  a  human  analyst  to  resolve  ambiguous  contact-to- 
track  associations  (the  ANALYST  module). 

•  updating  and  managing  a  database  of  the  status  of  all  contact  reports  and  tracks 
(the  Locational  Data  Base  Manager  and  MTST  modules). 


The  data  flow  through  these  modules  is  shown  in  Figure  1.   A  detection  is  made 
by  a  sensor  and  a  report  is  initiated.  After  the  applicable  delays  have  elapsed,  the 
contact  report  travels  along  the  communications  network  to  the  fusion  center, 
experiencing  additional  delay  at  each  node. 

Upon  arrival  at  the  fusion  center  the  INPUT  module  enters  it  into  the  Locational 
Data  Base  Manager  (LDBM)  module  and  passes  it  to  the  ASSOC  module.   There  the 
contact  is  compared  to  the  existing  tracks  on  the  basis  of  geofeasibility.  Those  tracks 
with  which  the  contact  could  be  associated  are  evaluated  by  a  statistical  comparison  of 
a  measure  of  correlation  (MOC)  with  a  preset  threshold.   If  the  contact  fails  to  trigger 
an  association  with  any  track  it  is  passed  to  the  START  module.  If  the  contact 
triggers  a  single  association,  or  if  one  track  MOC  exceeds  all  others  by  a  preset 
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Figure  2:  ASSET  contact  data  flow. 

margin,  it  is  unambiguously  associated  with  that  track  [Ref.  4:p.  3].  Ambiguous 
association  is  resolved  by  the  ANALYST  module,  which  has  a  given  probability  of 
making  the  correct  association  (.7)  and  an  associated  exponentially  distributed  delay 
time  spent  making  the  decision  (mean  of  .2  hours)  [Ref.  3:p.  11;  Ref.5: ACT- 
ANALYST]. 
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The  contacts  which  are  unresolved  are  passed  to  START  which  checks  the  value 
of  the  flag  single-report-to-track.  If  the  value  is  true,  START  initiates  a  new  track 
based  on  the  contact,  if  it  is  false,  the  contact  is  passed  to  the  multiple  contact  track 
initiation  module  CLUSTER.   This  represents  a  simplification  of  the  actual  OBU- 
ACT  START  module  which  uses  comparison  to  a  set  of  criteria  and  human  analyst 
interaction  to  determine  if  a  new  track  is  to  initiated  [Ref.  6:p.  3-4].   Within  the 
context  of  the  simulation,  this  simplification  makes  sense  as  two  classes  of  detecting 
sensors  are  represented;  those  that  report  continuous  contact  observations  such  as 
SOSUS  or  other  fixed  area  sensors,  which  are  processed  at  the  fusion  center;  and 
those  which  perform  platform-level  analysis  and  periodically  report  the  tracks  they 
hold,   such  as  submarines.   Platforms  in  the  first  category,  unless  possessing 
unusually  low  false  alarm  rates,  would  not  initiate  a  track  based  on  a  single  report, 
while  those  in  the  second,  would. 

Unassociated  contacts  which  remain  after  passing  through  the  ASSOC  and 
START  modules  are  processed  by  the  CLUSTER  module.   Here  subsets  of  all 
unassociated  contacts  are  evaluated  for  geofeasibility  and  constant  course  and  speed 
likelihood.  If  the  value  of  the  constant  course  and  speed  likelihood  exceeds  a 
threshold  value  then  the  subset  of  unassociated  contacts  is  used  to  initiate  a  new  track 
[Ref.  7:p.  1-3].   Contacts  which  are  not  associated  after  passing  through  ASSOC  and 
CLUSTER  carry  counters  which  increment  each  time  they  pass  through  the  two 
modules.   When  these  counts  exceed  maximum  limits  (hardwired  at  7  in  CLUSTER 
and  at  5  in  ASSOC)  they  are  dropped  from  the  LDBM  [Ref.  5:  OBU-ACT]. 
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When  a  new  track  is  initiated  or  a  contact  is  associated  with  a  track,  the  MTST 
module  is  called  to  filter  the  track  data  and  produce  an  optimal  estimate  of  the  target's 
location  and  velocity  along  with  the  error  covariance  matrix  describing  the  quality  of 
the  estimate.   The  full  capability  tracker  in  MTST  contains  the  matrix  versions  of  the 
Kalman  filter  equations.   These  equations  determine  the  mean  target  position  at  a 
future  point  in  time  and  the  covariance  matrix  representing  the  error  in  that 
prediction.   The  elements  of  this  covariance  matrix  can  be  recast  as  ellipses 
representing  two  a  or  86  percent  containment  regions  for  the  target's  possible  position 
and  velocity.   A  complete  description  of  MTST  is  the  subject  of  Chapter  n. 

The  foregoing  process  of  correlating  contacts  to  tracks  and  filtering  those  tracks 
to  provide  an  estimate  of  target  position  occurs  at  each  instance  of  a  data  fusion 
center  object.  When  formulating  an  architecture,  the  fusion  center's  role  must  be 
carefully  designed.   Since  no  contact  reports  flow  out  of  a  fusion  center,  it  is  not 
possible  to  directly  fuse  the  pictures  at  two  or  more  fusion  sites  into  a  higher  echelon 
center.   The  individual  sensors  originating  the  detection  reports  must  send  duplicate 
reports  to  all  fusion  centers  which  will  be  using  the  report  in  determining  a  tactical 
picture.   The  fusion  center  also  keeps  separate  instances  of  the  tracker  to  process 
surface  contacts  and  subsurface  contacts.  Together  with  the  START  processing  of 
single  contact  to  track  decisions,  this  assumes  that  a  great  deal  of  contact  processing 
is  going  on  below  the  level  of  the  fusion  center. 

Assuming  that  reports  arriving  at  the  fusion  center  are  pre-processed  justifies 
several  assumptions  that  the  DETECTION-REPORTER  object  makes  about  the 
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reports  it  generates.   Most  important  to  streamlining  filter  operation  is  the  assumption 
that  all  submarine  probability  areas  (SPA's)  are  circular.   While  this  represents  a 
considerable  inaccuracy  for  single  contact  line-of-bearing  detections,  if  one  assumes 
the  contact  report  represents  the  accumulation  of  a  significant  number  detections,  the 
circular  error  assumption  is  more  reasonable.   Several  types  of  detections  are  resolved 
in  ASSET  without  explicitly  modelling  the  geometry  of  the  searcher  and  target  at  the 
time  of  detection  (MPA  and  Fixed  Area  Sensors  specifically).   The  calculation  of  the 
size  and  orientation  of  explicit  AOU  ellipses  where  possible  would  slow  execution  and 
add  considerable  complexity  to  the  program,  requiring  even  more  specific  user  input 
sensor  parameters. 

The  assumption  that  this  circular  AOU  size  is  independent  of  range  is  more 
difficult  to  defend,  in  light  of  an  expected  linear  relationship  between  the  standard 
deviation  of  the  position  error  and  range.   That  assumption  on  ASSET'S  part  can  be 
remedied  easily,  for  the  cases  where  range  is  explicitly  determined,  by  expanding  the 
radius  of  a  sensor's  AOU  linearly  with  range  based  on  a  user  input  bearing  error. 
This  would  make  the  process  of  tabulating  preprocessed  variance  data  impossible, 
however,  as  the  sensor  AOUs  would  no  longer  be  constant.   A  good  compromise 
involves  defining  three  AOU  sizes  for  contacts  at  short,  medium  and  convergence 
zone  ranges.   This  would  allow  for  some  representation  of  the  range  effect  on  AOU 
size,  with  acceptable  complication  of  the  variance  pre-calculation  procedure. 

The  other  tracking  assumption  ASSET  makes  is  that  the  covariance  of  the 
velocity  can  be  represented  by  the  steady-state  stationary  cross-covariance  of  the  IOU 
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process,  independent  of  the  sensor  characteristics.   This  assumption  is  used  to 
initialize  the  covariance  matrix  of  a  new  track.   This  method  is  correct  for  a  position- 
only  contact  which  gives  no  information  about  velocity.   Given  position  and  velocity 
information  however,  this  is  only  accurate  if  the  sensor  is  so  inaccurate  that  the  IOU 
process  velocity  variance  is  small  compared  to  the  sensor  velocity  variance  [Ref.  8:p. 
5-20]. 

This  is  not  true  in  most  cases  however,  with  the  result  being  that  the  filter  re- 
initializes the  velocity  update  each  time  a  contact  is  processed,  never  including  any 
representation  of  the  velocity  accuracy  of  the  sensor  involved.  Thus,  the  velocity 
variance  is  based  solely  on  the  IOU  model  and  not  on  the  variance  of  the  velocity 
reported  by  the  sensors.   This  assumption  makes  the  tracker  more  responsive  to 
course  changes,  but  it  also  artificially  inflates  the  SPA  size  of  accurate  contact 
streams  because  it  ignores  the  sensors  velocity  error,  which  may  be  significantly  less 
than  the  IOU  value. 

Given  these  assumptions  and  the  noise  characteristics  of  the  IOU  process,  the 
filter  implementation  can  be  streamlined  to  better  match  the  data  provided.   This 
process  is  detailed  in  the  following  chapters  which  will  describe  the  IOU  process  and 
the  formulation  of  the  general  and  ASSET-specific  filter  equations. 
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n.     THE  MANEUVERING  TARGET  STATISTICAL  TRACKER 

Maneuvering  Target  Statistical  Tracker  (MTST)  is  the  name  given  to  a  class  of 
Kalman  filter-based  estimation  and  prediction  algorithms  utilizing  the  Integrated 
Ornstein-Uhlenbeck  process  to  model  the  statistics  of  target  motion.  To  successfully 
track  a  maneuvering  target,  the  details  of  the  target  path  should  be  statistically 
predictable.   How  accurately  an  automated  tracker  like  MTST  performs  its  task  is 
directly  related  to  how  well  the  tracker's  target  motion  model  describes  the  targets 
actual  random  motion. 

The  random  tour  model  [Ref.  9]  can  be  used  to  describe  the  statistical  properties 
of  a  target  moving  at  constant  speed  which  makes  random  course  changes  at  times 
separated  by  intervals  chosen  from  an  exponential  distribution.   This  provides  a 
reasonable  estimate  of  the  type  of  motion  expected  of  a  patrolling  submarine  target. 
Unfortunately,  the  target  distribution  generated  by  a  random  tour  is  not  normal  and 
cannot  be  used  directly  with  a  Kalman  filter  [Ref.  8:p.  2-10].   A  process  is  required 
which  produces  normal  distributions  which  best  approximate  the  mean  and  variance  of 
the  random  tour.   This  is  the  IOU  process.   The  mathematical  development  which 
follows  is  summarized  from  the  more  complete  treatment  in  Reference  8,  the  results 
of  which  are  in  agreement  with  References  5,  10  and  13. 
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A.      THE  INTEGRATED  ORNSTEIN-UHLENBECK  PROCESS 

The  Integrated  Ornstein-Uhlenbeck  Process  (IOU)  is  based  on  a  first  order 
stochastic  linear  differential  equation  which  approximates  the  higher  order  non-linear 
deferential  equation  which  defines  the  motion  of  a  randomly  maneuvering  target.   It  is 
a  member  of  a  class  of  equations  known  as  Langevin  equations  and  has  the  form: 

—u{t)  =  -pu(t)dt  +  adw(t).  (2.1) 

at 

The  first  term  represents  a  deceleration  caused  by  a  resistive  force  proportional  to  the 
velocity  u(t).  The  random  forces  acting  on  the  particle  are  represented  by  the  second 
term,  where  w(t)  is  a  Gaussian  white  noise  process. 

The  stochastic  nature  of  this  equation  makes  it  possible  to  find  only  the  statistics 
of  the  distribution  of  the  solutions,  rather  than  the  specific  solution  itself.   A  Gaussian 
w(t)  produces  a  Gaussian  velocity  process  u(t)  which,  like  any  normal  process,  is 
defined  completely  by  its  first  and  second  order  moments.  The  result  is  the  velocity 
distribution  of  a  particle  which  is  undergoing  random  motion  similar  to  Brownian 
motion,  experiencing  random  instantaneous  accelerations,  but  whose  velocity  is 
damped  by  a  spring-like  effect  which  constantly  accelerates  the  particle  in  the 
direction  opposite  its  velocity  at  a  rate  proportional  to  that  velocity.   The  position 
variance  of  such  a  particle  is  unbounded  over  time,  but  the  velocity  variance  is 
limited  by  the  damping  coefficient  0.  The  result  is  a  velocity  distribution  that  is 
normal  with  a  limiting  variance  given  by: 
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limV{n(0}  =  hmE{[u(t)-E{u(t)}]2}  =  ~,  P-V 

and  a  mean  for  any  given  /  given  by: 

E{u(t)}  =  e^E{u(0)}.  0-3) 

The  mean  and  variance  expressions  for  the  position  of  an  object  whose  motion 
is  governed  by  an  IOU  process  completely  specify  its  positional  distribution  over 
time.  In  order  to  show  that  the  IOU  process  approximates  the  same  motion  as  a 
random  tour,  the  radial  position  distributions  for  the  two  stochastic  processes  will  be 
shown  to  have  the  same  variance. 

B.   EQUIVALENCE  OF  THE  RANDOM  TOUR  AND  IOU  PROCESSES 

The  variance  of  the  radial  distance  from  (x(0),y(0))  is: 

E{R2(t)}  =  E{[x(t)-x(0)]2+\y(t)-y(0)?}.  a'4) 

For  the  Random  Tour,  the  following  holds: 

E{R\t)}  =  ^[at-Ue-<"]f  (2.5) 

or 

where  V  is  the  speed  of  the  randomly  touring  particle  and  a  is  the  mean  number  of 
course  changes  per  unit  time.  The  corresponding  result  for  the  IOU  process  is: 
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£{**(,)}  =  2£{fit-l+e'"],  (2.6) 

P 

where  a  is  the  scale  factor  for  the  random  acceleration  and  0  is  the  damping  constant 

on  the  velocity,  as  in  Equation  (2.1). 

These  two  equations  can  be  made  identical  given  the  following  relationships 

hold: 

0=a  f£  =  V2  (2.7) 

Thus,  a  Random  Tour  with  parameters  a  and  V  can  be  approximated  by  an  IOU 
process  with  parameters  |8  =  a  and  a  =  Va.m.   It  is  worth  noting  that  since  a  normal 
distribution  is  completely  specified  by  its  mean  and  variance,  the  IOU  process 
represents  the  best  normal  approximation  to  the  Random  Tour  [Ref.  10]. 

In  ASSET  the  IOU  parameters,  0  and  a  are  hardwired  to  represent  a  target 
conducting  a  random  tour  at  ten  knots  with  an  average  time  between  course  changes 
of  four  hours.   These  parameters  should  match  the  actual  motion  of  the  platforms 
modelled  in  a  given  scenario.   Since  the  user  can  choose  whatever  target  motion 
parameters  he  wants  for  each  individual  platform,  the  parameters  embedded  in 
ASSET  may  disagree  considerably  with  the  actual  target  motion.  This  disagreement 
causes  the  IOU  prediction  to  consistently  lead  or  lag  the  target's  mean  position, 
depending  on  whether  the  model's  speed  and  course  change  rate  correspond  to  a 
velocity  faster  or  slower  than  the  target's  speed,  respectively.   This  causes  the 
position  of  the  center  of  the  resulting  SPA  to  be  offset  from  the  true  target  position. 
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Even  when  the  tracker  parameters  exactly  match  the  target  motion,  a  considerable  lag 
occurs  if  the  time  between  measurements  is  not  significantly  greater  than  the  time 
between  maneuvers. 

C.      OTHER  MOTION  MODEL  OPTIONS. 

The  IOU  model  used  by  MTST  damps  out  the  mean  speed  of  the  target 
exponentially  over  time  to  account  for  target  maneuvers.   No  matter  how  small  the 
time  interval,  the  process  reduces  the  mean  speed  by  the  proper  amount.   If  the 
interval  between  detections  is  consistently  shorter  than  the  maneuver  interval,  the 
damped  velocity  is  used  to  predict  the  position  of  the  constant  velocity  target  and  the 
tracker  develops  a  lag.   One  remedy  for  these  inconsistencies  would  involve 
increasing  the  complexity  of  the  tracker  by  using  adaptive  methods  to  more  closely 
match  the  model  parameters  to  a  given  target's  track.   An  adaptive  filter  recognizes 
changes  in  the  targets  motion  and  compensates  for  it  by  changing  the  parameters  of 
the  motion  model. 

Desiring  to  reduce  the  filter's  complexity  however,  instead  of  explicitly 
modelling  an  adaptive  tracker,  the  use  of  such  a  tracker  can  be  approximated  using 
the  IOU  process  as  a  basis  for  the  size  of  the  SPA,  but  assuming  the  model  has  some 
adaptive  properties  to  position  it  more  accurately.   This  can  be  accomplished  using 
the  existing  IOU  model  to  estimate  the  SPA  variance  and  by  assuming  the  tracker's 
position  prediction's  are  distributed  about  the  target's  true  position  in  a  circular 
normal  fashion.   This  is  an  optimistic  assumption,  but  captures  the  fact  that  real  world 
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trackers,  using  more  complete  contact  data  than  ASSET  and  an  adaptive  motion 
model,  can  achieve  more  accurate  position  and  velocity  estimates.   Examinations  of 
such  trackers  can  be  found  in  References  14,  15  and  16.   Reference  14  in  particular 
demonstrates  a  passive  tracker  that  delivers  outstanding  accuracy  while  the  contact  is 
on  a  constant  course  and  speed,  but  has  difficulty  maintaining  track  during  difficult 
target  maneuvers. 

The  benefits  of  adaptive  tracking  are  not  applied  to  the  SPA  size  in  this  case,  as 
the  ASSET  tracker  is  not  an  adaptive  one,  but  are  assumed  in  simplifying  how  its 
position  is  determined.   This  technique  introduces  a  moderate  degree  of  error  in  SPA 
size  by  not  precisely  simulating  the  mechanics  of  distributing  the  position  of  the 
center  of  the  SPA.   It  also  does  not  take  into  account  the  filter's  velocity  components. 
Reference  15  provides  methods  for  countering  these  effects  through  noise  adaptation 
and  correlated  maneuver  gating  while  Reference  16  uses  bias-sensitive  maneuver 
detection  and  Kalman  gain  adaptation.    Another  option  used  to  improve  the  IOU 
model  is  the  use  of  a  dual  velocity  system  which  combines  a  short  term  IOU  process 
combined  with  a  long  term  one.   This  type  of  model  can  reproduce  a  variety  of 
motions  depending  on  the  weighting  factors  given  to  the  short  and  long  term 
components.  The  resultant  Kalman  filter  has  six  states  and  a  6  x  6  co variance  matrix. 
This  is  the  implementation  of  the  IOU  process  currently  used  in  the  TOMAHAWK 
fire  control  system  for  surface  ship  motion  modelling  [Ref.  8].  These  are  but  a  few 
of  the  schemes  more  complex  trackers  can  use  to  improve  motion  modelling.   The 
single  velocity  IOU  process,  while  not  adaptive,  comes  close  to  satisfying  the 
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necessary  conditions  and  since  it  is  the  basis  of  the  ASSET  tracker  will  be  used  for 
comparison. 

D.      THE  KALMAN  FILTER 

In  order  to  estimate  target  position,  a  method  of  combining  the  contact 
position  and  its  covariance  with  a  prediction  of  that  position  and  covariance  based  on 
the  IOU  process  is  necessary.   The  Kalman  Filter  provides  such  an  estimate,  provided 
the  process  model  is  linear  and  time  invariant,  which  the  IOU  process  is.    The 
Kalman  filter  is  derived  from  the  least-squares  estimation  theory  of  Gauss  and 
Legendre  [Ref.  ll:p.  7].   This  technique  sought  the  most  probable  value  of  an 
unknown  quantity  based  on  a  series  of  measurements  containing  unknown  errors. 
More  formally,  finding  the  most  probable  estimate  of  x  based  on  the  matrix  equation: 

Z„  =  H?  *  V,  (2.8) 

where  vector  zk  represents  a  measurement  taken  at  a  discrete  time  k,  based  on  an 
observation  matrix  H  revealing  information  about  one  or  more  state  parameters  of  the 
true  state  vector  x  but  corrupted  by  a  noise  vector  v.  This  technique  was  generalized 
by  Kalman  to  apply  to  linear  filter  theory,  producing  a  filter  which  produced  a  best 
linear  estimate  of  x  [Ref.  12]. 

The  resulting  filter  is  based  on  two  mechanisms  operating  sequentially,  one 
predicting  the  state  and  covariance  for  a  future  time  and  the  other  combining  this 
prediction  with  a  noisy  measurement  of  the  target's  state  taken  at  that  time.  This 
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should,  given  correct  motion  modelling  and  noise  parameters,  produce  a  state  and 
covariance  which  more  accurately  reflects  the  true  target  state  than  either  the 
prediction  or  the  measurement  alone. 

These  mechanisms  are  represented  by  two  systems  of  matrix  equations.  The 
variables  involved  are  defined  in  Table  1.   The  prediction  equations  for  the  state  x 
and  covariance  E  are: 

state  prediction:    it+1  =  <l>xk+nw,  (2.9) 

covariance  prediction:    Et+1  =  4>Lk<l>T  +  Qk,  (2.10) 

where  <f>  represents  the  transition  matrix  from  the  dynamic  system  model,  Q 
represents  the  variance  inherent  in  the  dynamic  model  and  ^  represents  the  mean  of 
the  noise  term  w  from  the  dynamic  model,  which  is  assumed  to  be  zero  in  ASSET. 
The  equations  for  combining  the  predicted  and  measured  states  and  variances  require 
the  computation  of  a  Kalman  gain  matrix  which  provides  the  weighting  factors  for  the 
combination.  These  equations  which  update  the  filter's  estimates  are: 

Kalman  Gain:    KM  =  tk.lHT[Htk.lHT  +  J^J"' ',  (2J1) 

(2.12) 
state  update:    xM  =  xM  +  Kk.x  [zM  -  fiv  -  HZM\ 

covariance  update:    Et+,  =  \I-KMti\tMi  (2.13) 

where  R  represents  the  variance  of  the  noise  inherent  in  the  measurements,  Z  is  the 
measurement  itself,  y^  is  the  mean  of  the  noise  inherent  in  the  measurement  (also 
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Table  I:  Kalman  Filter  Definition  of  Terms 


Kalman  Term: 

General  Form 

ASSET  Form 

State  vector  at  time  k: 

xk  n  x  1 

mean  4x1 

State  transition  matrix, 

time  k  to  time  k+1: 

<j>  n  x  n 

phi  4x4 

State  process  noise  mean: 

fMwnx  1 

assumed  zero 

State  error  covariance, 

at  time  k: 

^nxn 

var  4x4 

State  process  error 

covariance  at  time  k: 

jQknxn 

f  4x4 

Observation  matrix: 

H  m  x  n 

I  4x4 

Observation  at  time  k: 

Zk  m  x  1 

state  4x1 

Observation  error  mean: 

Hy  m  x  1 

assumed  zero 

Observation  error  covariance 

at  time: 

J?k  m  x  n 

covariance  4x4 

assumed  to  be  zero  in  ASSET)  and  /  is  the  identity  matrix  of  proper  size.   The  state 
and  covariance  updates,  assuming  H  (as  ASSET  does)  is  the  appropriate  identity 
matrix  can  be  written  as: 


state  update:    *t+1  =  J:M(tM'%t  -^"'zj. 


a- 14) 


covariance  update:    Ettl  =  (fcw"! +^ttl~1)  ,  (2-15) 

where  the  inverse  of  R^+i  and  Ek+1  are  the  weighting  matrices.   Thus  the  inverse  of 
the  noise  error  variances  are  the  weighting  factors  by  which  the  prediction  and 
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measurement  are  linearly  combined.   This  is  the  form  of  the  Kalman  equations  which 
ASSET  actually  uses.   The  advantage  of  this  form  is  the  ready  understanding  of  the 
underlying  mechanism  of  weighted  linear  combination.  The  disadvantage  is  the 
computational  necessity  of  performing  three  matrix  inversions  to  calculate  the 
variance  propagation.   Regardless  of  the  form  of  the  equations,  the  <£,  Q  and  R 
matrices  must  be  developed  to  complete  the  computational  form  of  the  filter. 

1.     Development  of  the  4>  Matrix. 

The  matrix  <t>  represents  the  transition  matrix  which  governs  the  dynamics  of  the 
state  between  two  discrete  times.  The  following  development  summarizes  the 
derivation  of  Reference  11,  Section  5.2.   The  differential  equations  defining  the  IOU 
process  and  velocity  can  be  combined  to  form  a  single  matrix  differential  equation 
which  takes  the  following  form  in  the  single  coordinate  direction  x: 


~x(t) 

"o 

l" 

~x(t) 

+ 

"o" 

_ii(0_ 

0 

0 

_u(t)_ 

a 

w(t) 


more  compactly:    X{t)  =  F(t)X(t)  +  Gw(t) 


(2.16) 


(2-17) 


where  Xfrlrepresents  the  system  state  vector  and  the  dot  over  a  symbol  indicates 
differentiation  with  respect  to  time.   Given  that  X(t)  has  the  value  X(tJ  at  time  t0  then 
X(t)  at  any  future  time  t>t0  is  given  by: 
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X(t)  =  *(r,OX(0  +  j  4>{t,r)  Gw(r)dr 


(2.18) 


where  ^(If,^  represents  the  transition  matrix  between  the  states  at  time  t0  and  time  t. 
Provided  that  F(t)  is  constant  in  time  and  assuming  that  t0  is  zero,  the  matrix  <f>  can  be 
shown  [REF.  8:p.  5-3,5-5]  to  have  the  form: 


4>(t)  = 


1 


0 


la-*-*) 


-# 


(2.19) 


This  represents  the  transformation  projecting  a  state  into  the  future  based  on  the  IOU 
model.   The  velocity  mean  decays  to  a  new  value  which  is  simply  the  old  value 
multiplied  by  e*\   The  position  is  moved  a  distance  equal  to  the  time  1//3  multiplied 
by  the  change  in  velocity: 


*(')*♦,  °x(f)k+±(u(t)k-u(f)M), 


(2.20) 


where: 


«(')*♦,  =  «(0*e-*. 


The  form  this  takes  when  expanded  to  cover  both  coordinate  directions  is: 


(2.21) 
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Ht)  = 


1 

0 

¥'- 

■*) 

0 

0 

1 

0 

1(1 -£?-*) 

0       } 

0 

0 

e-* 

0 

0 

0 

0 

e-* 

(2.22) 


where  f  is  time  elapsed  from  the  last  measurement  to  the  observation  time  of  the 
current  measurement.   It  is  important  to  notice  that  the  transition  matrix  depends  only 
on  the  IOU  parameter  /S  and  not  on  a.   Thus  the  IOU  velocity  damping  coefficient 
and  the  elapsed  time  determine  the  dynamics  of  the  mean  state  and  variance 
transition.   Appendix  A  contains  graphs  of  the  values  </>(l  2),  which  will  be  called  <£2, 
and  0(2), which  will  be  called  <£3,  take  over  various  time  intervals.   The  noise 
magnitude  coefficient  a  is  involved  in  the  noise  terms  R  and  Q. 

2.     Development  of  the  Q  Matrix. 

The  Q  matrix  represents  the  noise  introduced  by  the  white  Gaussian  process 
w(t).  As  discussed  above,  as  the  elements  of  this  matrix  get  large  the  filter 
increasingly  relies  on  the  measurement  data  and  ignores  its  predictions.   As  the 
elements  get  small  the  filter  increasingly  ignores  the  new  contact  data  and  the  errors 
in  its  predictions  compound  until  it  diverges,  completely  losing  track  of  where  its 
target  is.     The  definition  of  this  matrix,  the  process  noise  error  covariance,  is: 


Q  =  e[{Z-E(Z)}{?-E(Z)Y 
which  for  the  random  process  operating  in  (2.26)  becomes: 


(2.23) 


28 


Q  =  E' 


(  4>(t,r)  G  w(t)  dr        f  <f>(t,s)  G  w(s)  ds 

*t  sit 


(2.24) 


Since  the  mean  of  the  white  Gaussian  processes  w(t)  and  w(s)  are  both  zero  and  the 
variables  of  integration  r  and  s  are  equal,  this  reduces  to  the  form: 


Q    =     j<Kt,T)GGT<t>T(t,T)dT. 


(2.25) 


Substituting  the  appropriate  matrices  G  and  <t>  from  (2.25)  and  (2.29)  respectively  and 
performing  the  indicated  matrix  multiplication  results  in: 


Q  = 


i 


°L(l-e'W-*)e-i*-* 


0 


_(l-e  -«'->)c -«'-*> 


#e-m-* 


dr 


(2.26) 


for  the  single  coordinate  direction  x.  Assuming  that  to=0  and  recognizing  that  q(l,2) 
is  equal  to  q(2,l)  results  in  three  expressions  to  integrate  to  get  the  final  value  of  the 
Q  matrix  in  one  direction: 


Q  = 


qll       ql2 
q21       q22 


where    ql2  =  qll 


(2.27) 


qll  =  —  \2t -  1(3  - 4e "*  +  e -*»)1 


f2.28J 
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ql2  =  q21  =  ^(l-2e^  +  e^') 
q22-*L(l-e-»<) 


(2.29) 


(2.30) 


When  projected  into  two  coordinate  directions  the  full  representation  of  the  Q  matrix 


is: 


Q  = 


(2.31) 


qll  0  ql2         0  " 

0  qll         0  ql2 

ql2  0  q22         0 

0  ql2         0  #22 

Each  element  in  this  matrix  is  directly  proportional  to  o2  and  decreases  with 
increasing  /S.  Time  has  an  inverse  exponential  relationship  causing  the  noise  to 
increase  the  longer  the  interval  between  measurements.   At  t  equal  to  zero  the  noise 
terms  are  equal  to  zero  and  as  t  goes  to  infinity  : 


qll^t     «U^>     422^J- 


(2.32) 


2P' 


As  mentioned  previously,  the  variance  in  position  is  unbounded,  but  the  cross- 
covariance  and  velocity  covariance  are  bounded.  These  values  become  important  in 
developing  the  measurement  noise  covariance  matrix  R.   Appendix  A  contains  graphs 
of  these  values  over  various  time  intervals. 
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3.     Development  of  the  R  Matrix  and  Filter  Initialization. 

The  matrix  R  represents  the  covariance  of  the  measurement  noise  and 
accompanies  each  state  measurement  reflecting  its  precision.   This  value  may  vary 
with  time  in  the  ASSET  tracker  as  contact  reports  flow  in  from  a  variety  of  sensors. 
The  inverse  of  this  matrix  is  used  as  the  weighting  function  for  the  contact  report,  the 
larger  the  elements  of  R  the  less  effect  the  contact  has  on  the  filtered  update.   The  R 
matrix  also  plays  a  role  in  initializing  the  filter. 

When  a  target  is  first  detected,  there  is  no  t-1  state  for  the  filter  to  work  from  in 
making  its  prediction.   In  general,  a  time  t0  R  matrix  is  preset  to  initiate  the  filter. 
This  matrix  would  consist  of  the  2  x  2  positional  error  matrix  which  describes  an 
estimate  of  the  elliptical  AOU  about  the  contact  position: 


4 


{A 2  cos20)  +  (B2  sin20)      (B2  -  A 2)  sin0  cos0 
{B2  -  A 2)  sin0  cos0     {A 2  sin20)  +  (B2  cos20) 


(2.33) 


For  a  position-only  contact,  the  two  lower  right  diagonal  entries  would  be  set  equal  to 
the  stationary  variance  of  the  IOU  process  (2.2).  This  represents  the  maximum 
variance  the  velocity  could  possibly  have  based  on  the  motion  model  and  since  no 
velocity  information  is  coming  in,  the  filter  is  relying  solely  on  the  motion  model, 
making  this  a  value  a  reasonable  choice: 
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EQ.  2.33 

u       u 

(2.34) 
0       0 

0       0 


If  the  contact  is  a  position-and-velocity  contact,  the  procedure  is  more  complex.   In 
like  fashion  to  the  initialization  of  the  upper  left  2x2  portion  with  the  variance  in  the 
sensor  (its  elliptical  AOU),  so  must  the  bottom  right  2  x  2  be  initialized  by  the 
velocity  variance  of  the  sensor.   Since  the  IOU  process  is  predicting  the  velocity 
distribution  also,  this  variance  of  the  measurement  must  be  combined  with  the 
variance  of  the  IOU  process.   The  proper  combination  can  be  shown  [Ref.  8:App.  F] 
to  be  the  same  as  the  alternate  form  of  the  Kalman  covariance  update  equation  (2.14). 

This  method  is  not  used  by  the  tracker  in  ASSET.   The  ASSET  tracker  simply 
uses  the  initial  contact's  IOU  velocity  variance  as  the  basis  for  a  new  track,  without 
updating  it.  This  original  contact  report  is  used  as  the  t-1  state  and  covariance  when 
the  next  contact  associated  with  the  track  comes  in.  The  contact  reports  themselves 
also  differ  from  the  construction  described  above.  The  positional  variance  of  the 
sensors  are  all  circles  in  ASSET  so  the  2  x  2  matrix  which  occupies  the  upper  left  of 
the  R  matrix  is  actually  the  2  x  2  in  the  upper  left  of: 
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(2.35) 


The  velocity  variance  in  the  lower  right  hand  2  x  2  is  correctly  set  for  a  position-only 
contact,  but  this  is  also  used  for  a  position-and-velocity  contact.   Thus  the  R  matrix 
which  ASSET  uses  never  incorporates  the  error  in  velocity  inherent  in  the  sensor. 
Using  these  values  essentially  reinitializes  the  filter  each  time  a  measurement  is  made. 
This  is  desirable  for  position-only  contacts  as  it  prevents  divergence  of  the  velocity 
components  of  the  state,  however  for  position  and  velocity  contacts,  this  can  cause 
poor  velocity  accuracy  as  the  tracker  improperly  weights  the  velocity  components  of 
the  contact  report.   The  velocity  information  of  different  sensors  are  given  equal 
weight  regardless  of  the  relative  size  of  the  errors  in  the  reports  they  make.  The 
actual  velocity  variance  of  the  sensor  should  be  computed  and  used  to  set  the  value  of 
the  filtered  velocity  using  the  form  of  (2.15). 
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m.   IMPLEMENTATION  OF  THE  TRACKER  IN  ASSET 

Having  developed  the  components  of  the  Kalman  filter,  the  construction  of  the 
tracker  implemented  in  ASSET  can  be  discussed.   The  following  is  based  on 
Reference  5,  the  ASSET  source  code  for  version  1.0.  The  ACT-MTST  module 
consists  of  objects  which  perform  the  filter  calculations  described  above,  as  well  as 
objects  which  maintain  the  spatial  relationships  between  the  positions  involved  on  a 
spherical  earth. 

There  are  several  constants  which  are  defined  for  use  in  all  the  module's 
objects.   The  first  is  dtor  which  is  used  for  degree  to  radian  conversion.   Where 
trigonometric  functions  are  indicated  below,  the  actual  code  uses  dtor  to  convert 
angles  from  degrees  to  radians,  but  in  the  interest  of  brevity,  this  will  be  excluded. 
Second  is  a,  which  is  equivalent  to  0  in  the  IOU  process,  and  is  set  to  .25.   Third  is 
a  which  is  the  noise  scale  factor  from  the  IOU  process  as  described  above  and  set  to 
V50.  This  is  equivalent  to  the  square  root  of  a  velocity  multiplied  by  the  Random 
Tour  velocity.  Thus  the  tracker  predicts  the  target  position  by  assuming  a  random 
tour  is  being  conducted  at  a  speed  of  ten  knots  with  a  mean  interval  between  course 
changes  of  four  hours.   This  cannot  be  changed  by  the  user. 

Inputs  to  ACT-MTST  are  read  from  the  data  in  ACT-LDBM,  the  contact 
database  manager.   If  a  new  track  is  to  be  initiated,  the  next  MTST  object,  start-new- 
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track  sets  the  track  head  and  track  tail  state  and  covariance  matrices  equal  to  the 
contact  mean  and  covariance  using  compute-mean-co variance.   The  use  of  two  sets 
of  data,  a  head  and  a  tail,  to  allow  incorporation  of  out-of-sequence  contact  reports 
will  be  discussed  below.   The  object  compute-mean-covariance  is  also  called  by 
ACT-ASSOC  as  a  basis  for  the  spatial  measure  of  correlation  used  to  determine 
contact  to  track  correlation. 

Compute-mean-covariance  creates  two  arrays,  a  four  by  four  called  var,  which 
is  the  R  matrix,  and  a  four  by  one  called  mean,  which  is  4: 


mean  = 


0 

0 
spd  cos(90 -cse) 
spd  sin  (90  -cse) 


(3.1) 


and  var,  which  has  the  same  form  as  (2.34).  The  vector  mean  is  centered  at  the 
latitude  and  longitude  indicated  in  the  contact  data  list  and  x  and  v  velocity  computed 
from  the  contact  course  and  speed.   This  variance  matrix  however,  as  discussed 
above,  is  calculated  as  if  the  contact  report  consists  of  an  ellipse  with  minor  axis  B, 
major  axis  A  and  orientation  angle  6.  All  the  spatial  data  fields  that  are  generated  by 
the  Detection  Reporter  and  Generic  Sensor  contain  a  single  positional  uncertainty  and 
an  orientation  angle  of  zero.   Modifying  (3.1)  to  take  advantage  of  this  fact  results  in 
(2.35)  where  r  is  equal  to  the  standard  deviation  of  the  sensor's  target  position 
estimates.   This  alternative  form  returns  matrix  element  values  identical  to  those 
computed  by  (3.1),  given  the  actual  form  of  the  contact  reports  generated. 
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If  a  contact  is  associated  with  an  existing  track  rather  than  used  to  start  a  new 
one,  the  update-track  object  is  called  by  ACT-LDBM  to  perform  the  Kalman  filter 
operations.   The  remaining  objects  in  ACT-MTST  are  all  components  of  the  update- 
track  object  with  the  exception  of  the  object  cov-to-ellipse,  which  is  called  by  various 
display  and  graphics  objects  to  determine  the  orientation,  semi-major  and  major  axes 
of  the  ellipse  representing  a  covariance  matrix.   Since,  as  shown  above,  all  the 
ellipses  are  actually  circles,  this  object  is  another  anachronism  of  the  mismatch 
between  the  capabilities  of  the  tracker  and  the  contact  reports  actually  input  to  it.   It 
can  be  simplified  in  a  manner  similar  to  compute-mean  covariance  above  to  return 
the  radius  of  a  circle  instead  of  a  major  axis,  minor  axis  and  orientation  of  a  ellipse. 

Focusing  once  again  on  the  update-track  object,  it  first  defines  the  spatial  data 
parameters  of  the  contact  and  the  track  the  contact  is  used  to  update.   ASSET  uses 
four  data  structures  to  represent  the  data  associated  with  contacts  and  tracks: 


obu-contact:  consisting  of  contact  id,  receipt  time,  track  association,  sensor, 
categorization  (track  association  flags  "pinnedp"  and  "lockedp",  and  number  of 
passes  through  ACT-CLUSTER  and  ACT- ASSOC),  spatial  data  (the  obu-data- 
field  described  below),  altitude  (surface  or  subsurface)  and  HFDFp  (an  optional 
parameter  associated  with  HFDF  contacts). 

obu-data-fleld:  type  (position-only  or  position-and-velocity),  observation  time, 
latitude,  longitude,  AOU  orientation,  AOU  major  axis,  AOU  minor  axis, 
course,  course  uncertainty,  speed  and  speed  uncertainty. 

obu-track:  track  id,  number  of  contacts  associated  to  the  track  (maximum  of 
five),  state  and  covariance  of  the  track  head  (most  recent  contact)  and  track  tail 
(oldest  online  contact),  head  and  tail  spatial  data  (each  an  obu-data-fleld),  head 
and  tail  contact  id,  and  altitude  (surface  or  subsurface). 
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•  obu-state-field:  reference  latitude,  reference  longitude,  state  array  (offset  from 
reference  position  in  miles  N/S,  offset  from  reference  position  in  miles  E/W, 
N/S  velocity,  E/W  velocity)  and  covariance  array  (four  by  four  array  per  (3.1)). 


The  three  AOU  variables  can  be  eliminated  from  obu-data-field  and  replaced  by 
position-uncertainty,  the  radius  of  the  actual  AOU  circle  reported.  This  will  not 
affect  the  tracker's  calculation  and  reduces  the  memory  taken  up  by  contact  reports. 

An  obu-track  consists  of  up  to  five  contact  reports,  the  most  recent  of  which  is 
the  track  head  and  the  oldest  is  the  track  tail.  The  variables  required  for  performing 
the  filter  calculations  are  read  from  the  data  structures  described  above,  update-track 
then  checks  the  time  since  the  last  update.   If  the  new  contact  observation  time  is  later 
than  the  track  head  update  time,  the  contact  is  used  to  update  the  track  head  using  a 
set  sequence  of  object  calls.   This  sequence  contains  the  Kalman  filter  procedures  as 
well  as  two  routines  which  maintain  the  latitude  and  longitude  relationships  on  a 
spherical  earth  and  proceeds  as  follows: 


ioumotion:  Performs  the  IOU  prediction  step  and  updates  the  N/S  and  E/W 
offset  distances  in  the  track  head  state  to  the  time  of  observation  of  the  contact. 
These  offsets  effectively  perform  the  prediction  in  a  plane  tangent  to  the  earth's 
surface  at  the  reference  latitude  and  longitude  point  of  the  track  head. 

change-tangent:  Finds  the  relationship  of  the  iou  predicted  position  and  the 
contact  point  in  a  plane  tangent  to  the  earth  at  the  contact's  reference  point  vice 
the  track  head's  reference  point  used  to  compute  the  prediction  in  ioumotion. 
These  matrices  are  stored  as  state-cov-old.   Since  the  covariance  is  circular, 
there  is  no  need  for  the  sequence  which  rotates  it. 

make-state:  Takes  the  contact  data  and  produces  z*  which  is  called  simply 
mean  and  its  covariance  R  which  is  called  var.   These  two  matrices  are  stored 
in  state-cov-ctc. 
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•  filter:  Takes  state-cov-old  and  state-cov-ctc  and  uses  Equations  2.14  and  2.15 
to  compute  xk+i  and  1^+,.   The  result  is  the  Kalman  filter  update  performed  in  a 
plane  tangent  to  the  earth  at  the  contact's  reference  point. 

•  center-tangent:  The  latitude  and  longitude  of  the  position  produced  in  filter  is 
projected  back  onto  the  earth  from  the  tangent  plane. 

•  update-track-head:  The  new  updated  state  and  covariance  are  made  the  new 
track  head.   The  contact  used  to  perform  the  update  provides  the  track-head 
contact  data. 

•  update-track-tail:  Sorts  the  contacts  making  up  the  track  in  time  order,  oldest 
to  most  recent  (1  to  n).  If  there  is  only  a  single  contact,  it  is  designated  the 
track  tail.  If  there  are  five  or  fewer  contacts  and  the  oldest  contact  is  already 
the  tail  it  remains  so;  otherwise  the  oldest  contact  is  made  the  new  tail.   If  there 
are  more  than  five  contacts,  it  takes  number  contacts,  subtracts  five  and  the 
contact  whose  number  equals  the  result  is  used  to  update  the  track  tail  using  the 
procedure  detailed  above.  The  contact  used  to  update  the  tail  provides  the  track 
tail  contact  data.    All  contacts  older  than  the  tail  are  then  pruned  from  the 
track. 


Considerable  computation  time  can  be  saved  by  simulating  a  portion  of  the  above 
process,  rather  than  actually  carrying  out  all  the  steps.   A  revised  update  procedure 
based  on  using  the  updated  variance  (E^)  as  an  approximation  to  the  distribution  of 
updated  states  (*k+1)  is  proposed.   Rather  than  produce  the  contact  position  latitude 
and  longitude  randomly  using  the  sensor  variance  as  a  distribution  and  processing  it 
using  the  sequence  above,  the  latitude  and  longitude  are  chosen  by  computing  the 
updated  variance,  which  can  be  done  prior  to  scenario  execution.  When  it  is  needed, 
the  distribution  is  centered  on  the  target's  true  position,  and  used  to  generate  the 
updated  states  directly  by  random  draw.  This  new  procedure  analyzed  in  detail  in 
Chapter  IV. 
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Time  delays  incurred  at  communications  nodes  enroute  to  the  fusion  center  may 
cause  the  observation  time  to  fall  between  the  track  head  update  time  and  the  track 
tail  update  time.   In  that  case  a  modified  procedure  is  followed  . 

•  If  the  time  of  observation  for  the  contact  is  later  than  the  tail  update  time  or  if 
the  track  does  not  yet  consist  of  five  contacts,  then  the  contact  is  appended  to 
the  track  and  the  contacts  are  sorted  by  time  of  observation.   If  the  contact  is 
older  than  the  tail  and  the  track  already  consists  of  five  contacts,  it  is  ignored 
and  the  procedure  stops. 

•  If  the  observation  time  is  later  than  or  equal  to  the  tail  update  time,  the  tail  is 
used  as  the  oldest  contact,  if  the  new  contact  is  older,  then  compute-mean-var 
is  used  to  make  the  new  contact  the  oldest. 

•  The  oldest  contact  is  then  updated  using  the  procedure  described  above  using  the 
next  contact  in  the  track,  in  time  sequence,  as  the  "contact".   The  result  of  this 
becomes  the  "old  contact"  and  the  procedure  is  repeated,  updating  each  contact 
in  the  track  in  sequence  until  the  head  is  updated.   The  resulting  "forward 
filtered"  track  closely  approximates  the  track  head  which  would  have  resulted  if 
the  time-late  contact  had  been  received  in  sequence. 

The  calculations  performed  in  ioumotion  and  filter,  all  assume  that  elliptical 

AOU's  are  being  produced  by  the  Detection  Reporter  and  General  Sensor  objects. 

The  object  make-state  can  be  simplified  in  the  manner  as  compute-mean-cov  as  it 

performs  the  same  operations.   The  reduced  forms  of  the  equations  in  those  objects 

are  computed  below. 

A.      REDUCED  COMPLEXITY  IOU  PREDICTION 

The  matrix  equations  for  calculating  the  IOU  predicted  state  and  covariance  are 
given  by  Equations  2.9  and  2.10.   Using  the  symbolic  processing  capability  of  the 
MathCad  3.0  (copyright  1991  MathSoft,  Cambridge  MA)  computer  software  package, 
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the  matrix  operations  were  evaluated.   The  large  number  of  zero  terms  and  high 
symmetry  of  the  4>  and  Q  matrices  resulted  in  simple  expressions  for  the  elements  of 
the  state  and  variance  predictions.   Portions  of  the  calc-mot-mat  object  which 
compute  f(0,0),  f(0,2),  f(2,2)  for  positive  time  as  well  as  </>(0,2)  and  <i>(2,2)  can  be 
used  to  produce  the  terms  of  the  following  equations  used  to  compute  new-var: 


Xk^=newmean= 


E/W  posk^ 
N/S  posk+l 
EIW  ve/t+1 

— 

A 

- 

N/S  velk.x 

.V 

yk+phi(0,2)  Vyk 

xk+phi(0,2)  Vxk 

Vykphi(2,2) 

Vxkphi(2,2) 


(3.2) 


where: 


and: 


0(0,2)  =  1(1  -e-*) 
P 


and        <f>(2,2)  =  e'*; 


(3.3) 


tk.t=newvar= 


0      VaM     0      VbM 

VbM    o     VtM     o 


(3.4) 


where: 


and  with: 


**♦! 


Vak+2  (pW(0,2))  Vbk+phi(0,2)2  Vct+/(0,0) 

(Vbk+phi(0,2)  Vck)phi(2,2W(0,2) 

phiQa?  Vck+f(2,2) 


(3.5) 
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The/ terms  given  here  are  the  ASSET  versions  of  (2.36),  (2.37)  and  (2.38). 
B.      REDUCED  COMPLEXITY  KALMAN  UPDATE 

In  a  manner  analogous  to  that  used  in  the  preceding  treatment  of  the  IOU 
update,  the  matrix  equations  for  the  Kalman  update  were  evaluated.   The  result  is: 


Kalman  Gain  = 


Ka 


k*\ 


0 

Kb2 

0 


0 

Kat*i 
0 

Kb2M 


M      o 
0       Kbl 


k*\ 


KCk*\ 


0 


0        Kc 


k*\ 


where: 


(3.7) 
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Kb2k..  =1  [W.—  1 , 
Using  these  gain  values  the  corresponding  updated  variance  terms  are: 


'K*\ 


W>1M  =  KblM  *f 
Vb2k.x  =  Kb2M  Ir2 


Vct*i   ~  KCk*l  2Q 


(3.8) 


(3.9) 


The  updated  state  is  written  in  terms  of  residuals  which  are  the  differences  between 
the  contact  states  and  predicted  states: 
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Residuals  = 


Ry 

zy-xy 

Rx 

zx-xx 

RVy 

zVy-iVy 

RVx 

zVx-XVx 

y 

X 

Vy 
Vx 


$+KaRy+KblRVy 
x+KaRx+KblRVx 
Vy+Kb2Ry+KcRVy 
Vx+Kb2Rx+KcRVx 


(3.10) 


(3.U) 


The  reduced  form  equations  for  the  state  are  included  here,  but  are  unnecessary  if  the 
estimation  of  the  state  is  to  be  used.   The  accuracy  of  the  estimation  technique  is 
analyzed  in  the  following  chapter. 
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IV.  MODIFICATIONS  TO  THE  ASSET  TRACKER 

The  tracker  currently  used  in  ASSET  makes  two  major  assumptions  relating  to 
detection  modelling: 

•  Sensor  AOU's  are  circles  vice  ellipses, 

•  The  velocity  variance  inherent  in  the  sensor  is  ignored  when  making  the  filtered 
variance  calculation  for  position  and  velocity  contacts. 

Figure  2  shows  the  geometry  of  a  typical  filter  operation: 


Contact  Position 

Randomly  Drawn  From 

Distribution  Centered  on 
Target  »  True  Poaltion 


Now  Filtered  Position 
and  SPA 


r?\ 


Predicted  Position 
and  AOU 


New  Filtered  Poaltion 
and  8PA 


Contact  Position 


VI. ^ Target's  True  position 

Previous  Filtered  Position 
amd  SPA 


Figure  2:  Geometry  of  Filter  Update  Operation. 

The  equations  which  exactly  determine  the  position  of  the  center  of  the  SPA  are  given 
in  (3.12).  Given  that  the  simulation  knows  the  true  target  location  and  making 
several  assumptions  detailed  below,  it  is  proposed  to  bypass  the  Kalman  state  update 
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step  and  estimate  one  possible  position  of  the  SPA  based  on  the  relationship  between 
the  SPA  center  and  the  true  target  position.   If  the  true  position  is  distributed 
somewhere  in  the  SPA,  then  from  the  point  of  view  of  the  target,  the  SPA  center  is 
likewise  distributed  about  its  position.   This  allows  an  estimate  of  a  SPA  center  to  be 
generated  by  drawing  a  random  variable  from  a  distribution  with  standard  deviation 
equal  to  the  square  root  of  the  filtered  variance  in  position  with  a  mean  centered  on 
the  target's  true  position.   This  uses  the  SPA  "in  reverse"  to  find  a  possible  position 
for  its  own  center  based  on  the  targets  true  position.   Provided  the  necessary 
assumptions  are  met,  this  will  provide  a  reasonable  approximation  to  the  filtered  track 
over  a  series  of  observations.   It  also  allows  the  variance  values  to  be  precomputed 
and  drawn  from  a  table  vice  calculated  when  needed.   These  three  assumptions  are: 


•  The  mean  predicted  position  and  velocity  must  be  very  close  to  the  target's  true 
position  and  velocity  so  the  AOU  of  the  predicted  position  is  concentric  with  the 
contact  AOU,  or  nearly  so.   This  assumption  relates  to  the  motion  model  in  the 
tracker  simulated. 

•  It  is  sufficient  to  model  the  behavior  of  detections  in  the  mean,  rather  than 
explicitly  computing  them  individually.  This  assumption  relates  to  the 
interaction  between  the  correlator  and  the  tracker. 

•  The  distribution  of  filtered  states  has  its  mean  at  the  target's  true  position  and  a 
variance  equal  to  the  steady-state  value  of  Va  (3.10). 


The  first  assumption  sets  a  requirement  for  the  accuracy  of  the  motion  model. 
If  the  predicted  position  is  consistently  offset  a  considerable  distance  from  the  target's 
true  position  and  the  Kalman  gain  is  small,  then  the  center  of  the  filtered  position 
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distribution  will  be  drawn  off  the  target's  true  position.  This  can  be  shown  by 
rearranging  the  terms  in  Equation  (3.12)  to  yield  for  a  single  dimension: 

y  =  y  (l-Ka)  +  KaZy  +  Kbl  {ZVy-Vy)  (4-V 

The  filtered  position  y  approaches  the  contact  position  Zy  as  Ka  goes  to  one  and 
moves  off  to  the  predicted  position  $  as  Ka  goes  to  zero.   This  case  must  be  avoided 
to  prevent  significant  inconsistencies  between  the  distribution  of  actual  filtered 
positions  and  the  estimates  generated.   Working  to  keep  this  relationship  is  the  effect 
shown  in  (2.5)  which  acts  to  push  the  position  toward  the  observed  position  as  the 
time  interval  increases  for  a  given  observation  AOU  size  and  the  IOU  prediction 
variance  gets  large  compared  to  it. 

This  effect  is  not  easily  identifiable  in  the  Kalman  gain  equations,  but  increasing 
the  time  interval  increases  the  size  of  the  predicted  AOU  and  increases  Ka  for  a  given 
observation  AOU  size.   This  counteracts  the  tendency  of  the  predicted  position  to 
recede  to  the  previous  filtered  position  as  the  time  interval  increases  and  keeps  the 
required  relationship  intact  for  cases  where  the  time  interval  between  observations  is 
large.   The  maximum  error  will  occur  when  the  two  AOU's  are  nearly  the  same  size 
and  the  value  of  Ka  is  approximately  .5-.65.  This  assumption  is  also  involved  in  the 
computation  of  the  filtered  velocity,  where  it  is  assumed  that  the  difference  between 
the  predicted  and  observed  position  multiplied  by  gain  Kb2  is  small. 

The  second  assumption  requires  acceptance  of  a  lack  of  consistency  between  the 
correlator  and  the  tracker  models.    The  ACT-ASSOC  module  uses  the  Kalman  filter 
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for  calculation  of  the  spatial  measure  of  correlation  it  uses  to  associate  contacts  to 
tracks.   The  effect  of  using  equally  representative,  but  different,  implied  contact 
locations  by  the  correlator  and  tracker  to  do  contact  association  and  track  updating  has 
not  been  thoroughly  examined.  The  fact  that  the  filtered  position  produced  by  this 
method  does  not  relate  to  the  contact  position  used  for  track  association  does  not 
appear  to  be  a  serious  problem,  given  the  Monte  Carlo  nature  of  the  simulation.  The 
draw  used  by  the  correlator  can  be  looked  at  as  determining  whether  the  proper  track 
assignment  is  made,  a  false  assignment  made  or  no  assignment  made.   If  the  contact 
is  drawn  from  the  proper  distribution,  the  probability  that  any  single  draw  results  in 
any  of  these  occurrences  is  the  same  whether  or  not  the  tracker  goes  on  to  use  the 
contact  to  update  the  track  chosen. 

The  third  assumption  is  again  linked  to  the  accuracy  of  the  motion  model.  The 
result  of  making  this  assumption  is  shown  in  Figure  3.   The  first  part  is  essentially 
the  effect  of  assumption  one  being  true.   The  second  part  is  required  because  the 
variance  in  distribution  of  the  actual  filtered  position  is  affected  by  the  speed  of  the 
actual  target.   For  target  speeds  less  than  the  modelled  speed,  the  actual  distribution  is 
tighter  than  the  estimated  value.   The  opposite  is  true  if  the  target  speed  is  faster  than 
the  model.   The  implication  of  this  assumption  is  that  the  target  is  exactly  performing 
the  motion  modelled  and  the  tracker  does  not  respond  as  the  actual  Kalman  filter  does 
to  deviations  from  the  expected  motion.   Since  the  range  of  velocities  the  submarine 
targets  may  have  is  small,  this  deviation  is  typically  less  than  fifteen  percent  of  the 
SPA  size. 
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Estimated  Distribution  of  SPA  Centers 
Predicted  AOU 

Sensor  AO 
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Equal  to  the  filtered  Variance 


Target's  True  Position 


Previous  Filtered  Position 


Figure  3:  Modified  Filter  Operation  Geometry. 

A.      PROPOSED  MODIFICATIONS 

Implementation  of  these  changes  in  the  code  is  a  part  of  a  separate  effort  to  port 
ASSET  to  the  Sun  UNIX  workstation.   A  listing  of  the  applicable  code  is  contained  in 
Appendix  B.   The  update-track  object  described  in  Chapter  III  ia  the  primary  object 
to  be  modified  to  implement  these  changes.   Its  function  can  be  split  between  those 
used  to  precalculate  the  variance  values  and  those  used  to  estimate  the  velocity  if 
precalculation  is  desired.   This  involves  computing  the  three  steady  state  variance 
values  for  each  sensor  AOU  size,  for  representative  values  of  the  contact  interarrival 
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time.   Other  objects  to  be  changed  are  mentioned  below  and  included  in  the  Appendix 
B  as  well.    A  summary  of  these  modifications  include: 


•  ioumotion:  The  three  element  values  of  the  f  matrix  (f(0,0),  f(0,2)  and  f(2,2)) 
and  the  two  <j>  matrix  values  (0(2,2)  and  0(0,2))  are  computed  as  before,  but  not 
set  in  matrix  form.  The  values  are  used  to  compute  the  three  element  values  of 
the  predicted  variance  using  Equation  (3.6).  The  three  computed  values  will  be 
used  to  compute  the  steady-state  filtered  variance  values  which  are  stored  in  a 
table  for  future  reference.   A  table  of  the  number  of  iterations  needed  to  reach 
steady  state  versus  interarrival  time  is  detailed  in  Appendix  C. 

•  change-tangent:  This  object  is  not  required. 

•  make-state:  The  random  numbers  used  to  determine  the  contact's  latitude  and 
longitude  (from  the  contactPositionDraw  function  called  by  the  Detection 
Reporter)  can  be  retained,  but  as  discussed  above  this  seems  unnecessary.   The 
two  values  of  the  R  matrix  are  all  that  need  be  calculated  by  this  object. 

•  filter:  The  values  computed  in  ioumotion  are  used  to  calculate  the  three 
elements  of  the  filtered  variance  using  Equation  3.10  for  representative  values 
of  the  contact  report  interarrival  time.  A  graph  of  the  steady  state  gain  Ka 
versus  interarrival  time  is  contained  in  Appendix  D.   The  Kalman  gains  are 
computed  in  order  to  find  those  values  and  will  be  used  below  to  set  the  range 
of  operation  of  the  filter. 


This  set  of  routines  comprise  the  pre-processing  module  and  can  be  separated 
from  the  rest  of  MTST  so  it  computes  the  necessary  data  at  the  start  of  the 
simulation.  This  requires  computing  the  three  variance  values  for  33  interarrival 
times,  giving  approximately  .03  between  succesive  table  values  from  .03  to  .97. 
This  requires  100  values  be  tabulated  for  each  AOU  defined.  Retrieving  values  from 
a  data  array  of  perhaps  1500  elements  appears  to  be  significantly  faster  than 
calculating  the  values  individually  each  time  the  filter  is  called. 
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The  following  routines  perform  the  necessary  computations  at  the  time  the  filter 
operation  is  performed: 


center-tangent:  This  object  is  modified  to  call  the  function 
c  on  tact  Posit  ionDraw  and,  using  the  filtered  variance  drawn  from  the  table 
computed  above,  as  positional  variance  choose  a  center  for  the  SPA.  The 
distances  involved  to  not  require  the  use  of  the  spherical  earth  routine  in 
contactPosit ionDraw.   The  new  latitude  and  longitude  can  be  calculated  using 
1°  of  latitude=60.1077nm  and  1°  of  longitude=60.1077*cos(latitude)nm.  The 
error  in  this  estimate  is  less  than  one  percent  out  to  ranges  of  600nm  at  non- 
polar  latitudes.  The  filtered  velocity  can  be  estimated  by  the  following 
equation: 


Vy,.i  =  Vyte-«{\-K2)+K3ZVyt  (42) 

Vx^  =  Vxte'al{\-K3yK3ZVxt 

update-track-head:  The  position  and  velocity  computed  in  center  tangent  and 
the  new  variance  values  are  used  as  before  to  update  the  track-head.   The  state 
consists  of  the  position  latitude  and  longitude  and  the  corresponding  velocity 
components,  while  the  covariance  consists  of  the  three  steady-state  values  of  the 
filtered  variance. 

update-track-tail:  This  is  modified  in  the  same  way  as  update-track-head  to 
account  for  the  new  form  of  the  state  and  covariance. 


The  result  is  a  close  approximation  to  the  variance  values  obtained  from  the 
full  implementation  and  a  position  that  represents  the  result  of  some  combination  of 
predicted  position  and  contact  position  that  was  not  specifically  identified  and  operated 
on.  The  distribution  of  estimated  filtered  states  closely  approximates  that  of  actually 
computed  states.   The  specifics  of  this  comparison  are  detailed  below. 
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B.   ACCURACY  OF  THE  MODIFIED  TRACKER 

The  assumptions  underlying  this  new  procedure  introduce  deviations  from  the 
full  implementations.   As  discussed  above,  the  error  in  the  positioning  of  the  SPA 
about  the  target's  true  position  error  is  maximum  for  Ka  approximately  .5-. 65.   The 
mean  value  of  this  error,  assuming  the  worst  case  of  a  non-maneuvering  target,  can 
be  calculated  using  Equation  (3.3)  to  find  distance  from  the  old  filtered  position  to  the 
new  observation  position.   Since  the  observations  are  centered  on  the  true  target 
position,  the  mean  of  this  value  can  be  calculated  using  the  true  target  position.   The 
results  are  shown  in  Appendix  E  for  favorable  and  unfavorable  conditions.   From 
these  results  it  can  be  seen  that  the  error  is  excessive  for  targets  at  high  speed  if  the 
time  between  detections  is  very  long.  This  quickly  drops  off  however,  and  for  more 
reasonable  time  intervals  of  less  than  one  hour,  the  relative  error,  even  for  high  speed 
targets,  is  negligible. 

Assumption  three  must  also  be  valid  for  the  approximation  to  be  reasonable. 
The  validity  of  the  first  part  of  the  assumption  is  shown  above.   This  requires  that  the 
variance  of  the  distribution  of  filtered  positions  be  equal  to  the  filtered  variance. 
Appendix  F  presents  a  comparison  of  filtered  position  distributions  and  filtered  AOU 
sizes  for  various  target  speeds,  computed  using  the  model  in  the  Appendix  G.   The 
data  in  Appendix  F  is  the  result  of  modelling  48  hours  worth  of  tracking  data  or  no 
more  than  300  observations.   The  simulated  target  conducted  a  random  tour  consisting 
of  ten  course  changes  executed  at  time  intervals  with  an  exponential  distribution  and  a 
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mean  of  4  hours.  The  target's  speed  was  either  3  knots,  10  knots  or  30  knots.  The 
standard  ASSET  IOU  parameters  /3=.25  and  o=lO0m  were  used.  The  results  of  30 
different  iterations  were  averaged  to  determine  the  mean  values. 

The  results  indicate  that  for  long  time  intervals  between  contacts,  and  speeds 
near  the  modeled  one,  the  results  are  fairly  good.   Excessive  deviations  are  exhibited 
when  the  target's  speed  is  far  in  excess  of  the  modeled  speed.   Even  in  the  case 
where  the  actual  target  speed  was  ten  knots,  the  damping  of  the  velocity  is  apparent  in 
the  filtered  mean  speed.  When  this  is  combined  with  the  offset  error  described  above, 
errors  of  up  to  sixteen  percent  can  result.   While  this  is  a  large  error,  the  average 
value  this  error  takes  through  the  range  of  common  time  intervals  between 
observations  (.  1  hours  to  1  hour)  and  possible  contact  speeds,  is  closer  to  five 
percent. 

One  source  of  this  error  can  be  seen  by  comparing  (3.10)  and  (3.12).   While  the 
single  Kalman  gain  Ka  is  used  to  compute  the  value  of  the  variance,  the  actual 
distribution  is  computed  using  both  Ka  and  Kbl.   For  the  range  of  AOU  sizes 
associated  with  submarines  and  MPA  assets  (at  worst  25nm)  leaving  the  SPA 
unadjusted  amounts  to  an  error  in  standard  deviation  of  the  estimated  distribution  of 
less  than  lnm.  The  error  is  typically  conservative,  increasing  the  SPA  size  by  the 
error,  and  is  largest  when  the  time  interval  between  contacts  is  small.   This  produces 
errors  which,  while  a  significant  percentage  of  the  AOU  size,  are  actually  quite  small. 
A  notable  exception  is  the  decrease  in  the  standard  deviation  in  situations  where  the 
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target  speed  is  high,  counteracting  some  of  the  IOU  process'  tendency  to 
underestimate  the  targets  speed  in  this  case. 

The  failure  of  the  estimate  to  account  for  the  target's  velocity  causes  the 
resultant  track  to  be  erratic.   This  effect  worsens  as  the  distance  traveled  between 
observations  gets  small  with  respect  to  the  AOU  size.   The  actual  filter  exhibits 
consistant  errors  in  this  case,  while  the  estimate  zig-zags  wildly  around  the  true  target 
track.   Since  the  linearity  of  the  track  is  not  critical  to  any  part  of  ASSET  except 
possibly  the  correlator,  the  effect  is  minimal.   For  any  single  point  the  estimate 
produces,  a  track  could  exist  which  includes  it.  The  random  draw  combines  elements 
from  all  these  hypothetical  tracks  into  a  single  collection  which  has  the  proper 
pointwise  relationship  to  the  real  track.   Thus,  while  specific  combinations  of 
observation  AOU  size  and  observation  interval  result  in  poor  performance,  others 
result  in  improved  performance. 

The  performance  of  the  modified  filter  is  reasonably  consistent  with  the  original 
filter  operation.   A  3-5  nm  error  in  a  25  nm  AOU  does  not  significantly  impact  the 
cueing  or  search  efforts  of  prosecuting  platforms  and  tracker  performance  is  actually 
improved  for  cases  where  the  target  is  transiting,  moving  at  high  speeds  without 
maneuvering.   The  errors  present  in  the  modified  tracker  appear  to  have  no  negative 
impact  on  the  operation  of  the  rest  of  the  simulation  and  while  introducing  moderate 
errors  in  tracking  maneuvering  targets,  balances  that  with  improved  performance 
against  non-maneuvering  ones.   Appendix  G  contains  examples  of  the  original  filter's 
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tracking  ability  and  that  of  the  estimate  for  a  variety  of  sensor  AOU  sizes  and  mean 
interarrival  times. 

C.      LIMITING  THE  RANGE  OF  FILTER  OPERATION 

An  interesting  effect  apparent  from  Appendix  F,  is  the  limiting  nature  on  the 
filtered  position  distribution  of  the  contact  AOU  on  the  one  hand  and  the  predicted 
AOU  on  the  other.   This  situation  stems  from  the  same  effect  that  moves  the  filtered 
position  from  the  observation  position  toward  the  predicted  position.   This  can  be 
seen  directly  by  looking  at  the  values  of  the  Kalman  gains  as  a  function  of  range  for  a 
given  time  interval  between  observations.   Examples  of  these  relationships  are  shown 
in  Appendix  H.   As  can  be  seen  from  these  graphs,  if  the  observation  AOU  radius  in 
nautical  miles  is  limited  to  being  less  than  150  times  the  time  interval  in  hours,  the 
Ka  value  remains  above  the  "flat"  portion  of  the  curves,  approaching  zero. 
Observations  which  arrive  at  the  tracker  such  that  the  interval  since  the  last  contact  is 
too  small  compared  to  the  AOU  size  will  be  filtered  using  primarily  the  prediction 
information.  The  effect  of  the  contact  is  largely  ignored.   This  is  the  filter's  way  of 
saying  that  the  information  value  of  the  contact  is  not  high  enough  to  be  considered  of 
value  in  computing  the  filtered  position. 

Such  contacts  can  be  ignored  as  there  is  little  difference  between  the  predicted 
position  and  variance  and  the  resultant  state  and  variance.   A  lower  limit  of  lnm  on 
the  observation  is  hardwired  into  the  make-state  and  compute-mean-cov  objects.   If 
the  time  between  contacts  is  too  short,  little  is  gained  unless  the  AOU  is  tiny.  The 
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lnm  limit  implies  that  the  minimum  time  interval  between  contacts  is  around  .01 
hours  to  allow  the  contact  data  to  be  assimilated  into  the  track.   If  the  contact  makes 
no  contribution  to  the  track  update,  there  is  no  reason  to  perform  the  filter  operation. 
The  effect  of  this  can  be  seen  in  the  data  in  Appendix  F,  where  the  predicted  and 
filtered  values  are  very  close  for  small  time  intervals  and  large  observation  AOUs. 
Thus  in  dense  contact  environments,  with  contacts  arriving  at  the  fusion  center  at 
intervals  of  less  than  30  minutes,  any  contact  based  on  an  AOU  larger  than  150  times 
the  time  since  the  last  contact  in  hours,  can  be  ignored  and  not  filtered.   This  value  of 
150  corresponds  to  a  Kalman  Gain  limit  of  approximately  .1.   The  exact  value  is 
rather  arbitrary  so  long  as  it  not  too  far  into  the  portion  of  the  curves  that  change 
quickly  with  variation  in  observed  AOU  size.   This  could  be  made  a  user  input  to 
limit  the  resolution  of  the  tracker  and  speed  execution  time  during  scenario 
development  and  testing,  and  then  reduced  to  provide  more  accurate  track  information 
when  the  data  is  actually  collected. 

In  similar  fashion  tracks  which  have  not  been  updated  in  a  long  time,  do  not 
contain  enough  information  to  affect  contacts  that  arrive  with  AOUs  less  than  5-10 
times  the  time  interval  since  the  last  observation.   If  the  value  of  Ka  is  very  close  to 
the  "flat"  portion  approaching  one,  the  predicted  value  is  ignored  in  favor  of  the 
contact  data.   This  allows  similar  saving  of  filtering  operations  in  sparse  track 
conditions  with  accurate  sensors.   If  the  time  interval  between  contacts  is  expected  to 
be  greater  than  one  hour  and  the  contact  AOUs  are  small,  using  the  raw  contact  to 
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update  the  track  is  a  reasonable  approximation  if  the  observation  AOU  size  is  smaller 
than  five  times  the  time  between  contacts  in  hours. 

When  these  two  "gating"  techniques  are  utilized,  substantial  numbers  of  filtering 
calls  can  be  eliminated.   This  can  be  particularly  true  in  scenarios  covering  large 
ocean  areas  with  many  possible  targets.   During  the  initial  stages  of  the  campaign,  as 
search  assets  slowly  get  cues  for  very  quiet  enemy  tactical  platforms,  detections  will 
not  be  made  until  the  enemy  is  at  close  range.   This  will  result  in  observations  with 
relatively  small  AOUs  coming  in  at  long  time  intervals.   As  cues  are  dispatched  and 
friendly  forces  converge  on  the  enemy,  even  a  modest  number  of  platforms  soon 
begin  producing  contact  streams  that  arrive  at  rapidly  decreasing  intervals.   Putting 
limits  on  when  the  filter  is  invoked  will  prevent  the  simulation  from  wasting  efforts 
computing  a  SPA  that  does  not  really  use  much  of  the  contact  information. 
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V.   CONCLUSION 

From  a  discussion  of  the  construction  and  use  of  the  Kalman  filter  based 
Maneuvering  Target  Statistical  Tracker  it  has  been  shown  that  the  tracker 
implementation  in  ASSET  is  not  currently  matched  well  to  the  types  of  inputs  it 
receives.   This  results  in  the  calculation  of  many  duplicate  values,  the  performance  of 
unnecessary  matrix  manipulations  and  the  calculation  of  values  in  situations  when 
little  useful  information  is  gained.   Improvements  to  the  tracker  aimed  at  decreasing 
the  computation  time  which  were  discussed  include: 


•  Use  of  the  equations  detailed  in  (3.3)  through  (3.12)  to  compute  the  actual 
filtered  variance  and  position. 

•  Filter  the  variance  of  the  contact  distribution  about  the  target  prior  to  conducting 
a  random  draw,  thus  producing  a  position  drawn  from  a  distribution 
approximating  that  of  the  actual  filtered  positions  variance.   Since  the 
observations  play  no  part  in  calculating  this  variance,  appropriate  values  for  the 
sensor  AOUs  can  be  precomputed  for  a  range  of  expected  time  intervals. 

•  Limit  the  situations  where  the  filter  is  actually  used  to  those  which  will  result  in 
a  meaningful  amount  of  information  being  extracted. 

•  Using  planar  estimations  of  the  latitude  and  longitude  computations  rather  than 
perform  the  conversions  on  a  "spherical  earth". 


When  taken  together  these  modifications  greatly  reduce  the  complexity  of  the 
calculations  required  to  update  a  track  after  an  observation.   They  can  also  be  put 
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under  user  control  so  the  level  of  fidelity  can  be  adjusted  based  on  the  requirements 
of  the  architecture  and  the  stage  in  development  of  the  architecture  the  analyst  is  at. 
Areas  for  future  work  in  this  area  include: 

•  Formulating  a  faster  way  of  estimating  the  steady  state  gains  for  use  in  a  table 
of  pre-calculated  filtered  variances. 

•  Analysis  of  the  filtered  position  distributions  of  other  types  of  trackers  to 
determine  if  a  better  method  of  approximation  can  be  determined,  one  which 
better  accounts  for  target  velocity. 

•  Formulating  a  simple,  recursive  estimate  of  the  Kalman  gains  for  on-the-fly 
calculation.   An  attempt  at  this  by  using  limiting  values  as  the  time  interval  got 
small  yielded  very  good  steady  state  results,  but  divergence  and  chaotic 
behavior  for  small  AOU  sizes  ( <  lOnm  radius)  precluded  its  use  iteratively. 

•  An  analysis  of  the  algorithms  in  the  correlator  to  determine  if  all  of  its 
algorithms  are  necessary,  or  if  a  simple  probability  of  correct/false/no 
correlation  could  be  substituted.  If  appropriate,  determine  what  those 
probabilities  should  be. 

•  Modelling  how  the  AOU  size  should  change  with  respect  to  range  and  bearing 
from  the  searcher.   A  strict  linear  bearing  error  vs  range  is  a  simple  way  of 
solving  the  problem  of  constant  sensor  AOU  size.   The  real  situation  is  not 
quite  so  simple.  This  requires  computation  of  the  beam  sizes  and  beam 
distribution  of  a  towed  array  and  the  effect  of  TMA  techniques  on  determining 
the  position  course  and  speed  of  the  target.  This  would  allow  for  proper 
calculation  and  filtering  of  the  sensor's  variance  in  velocity  measurement  as 
well. 

The  last  example  opens  up  a  plethora  of  possibilities  for  analysis  of  the  parameters 

required  for  input  to  ASSET.   A  compendium  of  data  on  everything  from 

communication  network  modelling,  to  aircraft  maintenance  data  analysis,  to  platform 

on  platform  engagement  simulation,  and  so  on.  Well  researched  values  for  the 

present  and  expected  future  interactions  of  platforms  and  their  supporting 
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infrastructures  need  to  be  produced  if  ASSET  and  simulations  like  it  are  to  be 
available  to  the  average  user. 

ASSET  has  the  potential  to  be  a  very  effective  analysis  tool.   A  great  deal  of 
work  must  be  expended  to  provide  the  necessary  data  to  allow  its  practical  use.   In  its 
present  form  it  requires  so  much  expert  knowledge  on  so  diverse  a  collection  of 
topics,  it  is  doubtful  one  person  could  perform  meaningful  analysis  without  an 
inordinate  amount  of  time  spent  in  research.   Continued  evolution  in  the  functionality 
of  the  simulation  will  provide  increased  applicability  and  accuracy  of  results. 
Increased  streamlining  of  the  computational  algorithms  will  be  necessary  as  this 
occurs  to  keep  it  running  on  a  desktop  workstation  in  a  reasonable  period  of  time. 
Striving  for  improvement  in  performance  however,  without  providing  a  ready  source 
of  data  will  most  probably  preclude  widespread  use  of  the  simulation. 
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APPENDIX  A:  CHARACTERIZATION  AND  ESTIMATION  OF  <t>  AND  Q. 


The  following  equations  represent  the  elemets  in  the  transition  matrix,  phi,  and  the  IOU  process  noise,  Q. 
Approximation  equations  which  can  be  used  to  more  quicly  calculate  the  values  are  also  given. 


<*>- 


1  0  02     0' 

0  1      0  4>2 

0  0  <t>3     0 

0  0     0  <t>3 


q1      °     <12     0  1    t    s    1,10. .5000      timet   =    (t).01 


Q- 


0  ql  0  q2 
q2  0  q3  0 
0     q2      0     q3 


tome. 


n    =    1  ,10..  300 
V    =    10        a 


.25 


V-Ja 


1 


Actual  value: 


Estimated  Value: 


<f>2t       -  •  f  1  -  expf    a  •  time, 
1        a  L  i  l 


1     . 


4>2et  -   -y  •  timet  •  [  -  2  +a  •  timet 


4>3,   -    exp[  -  a  •  timet 


tf>3e,  =    1  +t  •  a  •  timet  •  [  a  •  timet  -  2 


*  ■  i 


time. 


-  •  [  [  2  •  [  1  -  exp[  -  a  •  timet  1 1  -  .5  •  [  1  -  exp[  -  a  •  2  ■  time, 


qle,   -    <r-\  time. 


1     1 
3"4-o-t»nc 


a2  r 


q2,   =    —z>-ri  .5  -  expr    a-  time  J  M  .5  expf  -  a -2 -time,  ;  q2e,       a*-  =  -!  time.  !*•!  1     a -time, 


j2-l-r  time.]2-  1     a- 


a 


q3,    =    .5-1     exp     a -2   time,1 
1  at 


1  i 


q3et   =      la-  time,  • ;     1  +a  •  time, 


The  graphs  which  follow  show  the  behavior  of  these  functions  over  various  time  intervals.   The  variable  "time" 
is  equal  to  the  elapsed  time  from  the  last  contact's  time  of  observation  to  that  of  the  present  contact,  in  hours. 
The  estimated  values  are  computed  using  the  first  two  significant  terms  in  the  series  expansion  of  the  exponential 
terms.   These  estimates  are  very  accurate  for  interarrival  times  of  two  hours  or  less. 
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APPENDIX  B:  APPLICABLE  SOURCE  CODE 

This  appendix  contains  portions  of  the  ASSET  source  code  which  would  be  modified 
to  implement  the  changes  described  above.   The  remarks  indicate  the  appropriate 
portions  which  would  be  deleted  or  changed.   This  is  not  an  all  inclusive  listing,  but 
contains  the  portions  most  affected  by  the  modifications. 

ACT-MTST  Wed,  Apr  25,  1M0  1 


OBJECT  ACT-MTST 

ACT-KTST  Is  an  Object  LISP  implementation  of  KTST,  the 
Manuevering  Target  Statistical  Tracker. 


(setq  act-mtst  (kindof  nil)) 


(defobfun  (exist  act-mtst)  unit-list) 

(usual-exist  lnit-list) 

(nave  'dtor  (/  pi  180)) 

(have  'alpha  0.2S) 

(nave  'slgma  isqrt  (/  100  2)))       ;  10  knot  target 
) 


(defobfun  (create-objeet-llnks  act-mtst)  (osis) 
(have  ' ldbm  (ask  osis  ldbm) ) 

) 


(defobfun  (start-new-track  act-mtst)  (uc  track-id) 
(let*  ( (mean-cov-list  (compute-mean-covarlance  uc) ) 
(mean  (first  mean-cov-list) ) 
(var   (second  mean-cov-list) ) 
(alt   (obu -con tact- altitude  uc) ) 
(uc-spatlal-data  (obu-contact-spatial-data  uc) ) 
(lat   (obu-data-f leld-lat  uc-spatlal-data) ) 
(lng   (obu-data-f leld-lng  uc-spatlal-data) ) 
(ctc-list  (list  uc)) 
track  state-cov 
) 
(setq  state-cov  (make-obu-state-field 

: ref-iat  lat   :ref-lng    lng 
: state   mean  :covariance  var)) 
(return-from  start-new-track 

(make-obu-track  :id  track-id 

:number-of -contacts  1 
: head-state-covarlance  state-cov 
: head-contact -id  (obu-contact-id  uc) 
: head-spatial -data  uc-spatial-data 
:tail-state-covariance  state-cov 
: tail-contact-id  (obu-contact-id  uc) 
: tall-spatial-data  uc-spatlal-data 
: altitude  alt 
: contacts  ctc-list 
) 
) 
) 
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ACT-MTST  Wed,  Apr  25, 1990 


(defobfun  (compute-mean-covariance  act-mtst)  (uc) 

(let*  ( (uc-spatial-data  (obu-contact-spatial-data  uc) ) 
(type  (obu-data-f ield-type  uc-spatlal-data) ) 
(lat   (obu-data-field-lat  uc-spatial-data)) 
(lng   (obu-data-f ield-lng  uc-spatlal-data) ) 
(brg   (obu-data-f ield-brg  uc-spatial-data) ) 
(maj   (obu-data-field-maj  uc-spatial-data) ) 
(min   (obu-data-field-min  uc-spatial-data) ) 
(cse   (obu-data-field-cse  uc-spatial-data) ) 
(cse-unc  (obu-data-f ield-cse-unc  uc-spatial-data) ) 
(spd   (obu-data-field-spd  uc-spatial-data)) 
(spd-unc   (obu-data-field-spd-unc  uc-spatial-data)) 
(mean  (make-array  4      : initial -element  0)) 
(vax   (make-array  '(4  4)  : initial -element  0)) 
stationary  smaj  smin  cosb  sinb  cosb2  sinb2 

) 

(setq  stationary  (/  (•  sigma  sigma)  (*  2  alpha))) 

(setf  (aref  mean  2)  (•  spd  (cos  (•  dtor  (-  90  cse))))) 

(setf  (aref  mean  3)  (*  spd  (sin  (*  dtor  (-  90  cse))))) 

(setf  (aref  var  2  2)  stationary) 

(setf  (aref  var  3  3)  stationary) 

(setq  smaj  (max  maj  1)) 

(setq  smin  (max  min  1)) 

(setq  cosb  (cos  (*  dtor  brg))) 

(setq  sinb  (sin  (*  dtor  brg))) 

(setq  cosb2  (*  cosb  cosb)) 

(setq  sinb2  (*  sinb  sinb)) 

(setq  smaj2  (*  0.25  smaj  smaj)) 

(setq  smin2  (*  0.25  smin  smin)) 

(setf  (aref  var  0  0)  (+  (*  smin2  cosb2)  (*  smaj2  sinb2) ) ) 

(setf  (aref  var  0  1)  (*  (-  smaj2  smin2)  sinb  cosb)) 

(setf  (aref  var  1  0)  (aref  var  0  1)) 

(setf  (aref  var  1  1)  <♦  (*  smin2  slnb2)  (*  smaj2  cosb2))) 

(return-from  compute-mean-covariance  (list  mean  var) ) 


(defobfun  (update-track  act-mtst)  (uc  track) 

(let*  ((uc-spatlal-data  (obu-contact-spatial-data  uc) ) 
(id  (obu-contact-id         uc) ) 

(type  (obu-data-f leld-type  uc-spatial-data) ) 
(lat  (obu-data-field-lat  uc-spatlal-data) ) 
(lng  (obu-data-f leld-lng  uc-spatial-data) ) 
(brg  (obu-data-f ield-brg  uc-spatial-data) ) 
(maj  (obu-data-field-maj  uc-spatlal-data) ) 
(min  (obu-data-field-min  uc-spatlal-data) ) 
(cse  (obu-data-field-cse  uc-spatial-data) ) 
(cse-unc  (obu-data-fleld-cse-unc  uc-spatial-data) ) 
(spd  (obu-data-field-spd  uc-spatial-data)) 
(spd-unc  (obu-data-field-spd-unc  uc-spatial-data)) 
(obs-time   (obu-data-field-obs-tlme  uc-spatial-data) ) 
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(track-spatial -data    (obu-track-head-spatial-data  track)) 
(tail-spatial-data      (obu-track-tail-spatlal-data  track) ) 
(last-update-time    (obu-data-field-obs-time  track-spatial -data) ) 
(tail-update-time    (obu-data-field-obs-time  tail-spatial-data   ) ) 
(last-state-cov    (obu-track-head-state-covariance  track)) 
(tail-state-cov    (obu-track-tail-state-covariance  track)) 
(last-ref-lat    (obu-state-fleld-ref-lat   last-state-cov)) 
(last-ref-lng   (obu-state-fleld-ref-lng  last-state-cov)) 
state-cov-old  state-cov-ctc  state-cov-new  delta-time 
contacts  time-old  mean-var-list  pos  etc  ctc-sd  count 


) 


(setq  delta-time  (-  obs-time  last-update-tlme) ) 

(if  (>-  delta-time  0) 
(progn 

(setq  state-cov-old  (change-tangent  lat  lng 

(ioumotion  delta-time  last-state-cov) ) ) 
(setq  state-cov-ctc  (make-state  type  lat  lng  brg  maj  min  cse  cse-unc 

spd  spd-unc) ) 
(setq  state-cov-new  (center-tangent  (filter  state-cov-old 
state-cov-ctc) ) ) 

(setq  track  (update-track-head  track  uc  state-cov-new 
uc-spatlal-data  id)) 

(setq  track  (update-track-tail  track)) 
) 
(progn 

(if  (or  (>  obs-time  tail-update-time) 

(<  (obu-track-number-of-contacts  track) 
(ask  ldbm  contacts-from-head-to-tail) ) ) 
(progn 

(setq  contacts  (append  (obu-track-contacts  track)  (list  uc))) 
(setq  contacts  (sort  contacts  *' (lambda  (x  y) 

(<  (obu-data-field-obs-time  (obu-contact-spatial-data  x) ) 
(obu-data-field-obs-time  (obu-contact-spatial-data 


y) ))) )) 


(if  (>-  obs-time  tail-update-time) 
(progn 

(setq  state-cov-old  tail-state-cov  ) 
(setq  time-old      tail-update-time) 
) 


contacts) ) ) 

(first  contacts) ) ) ) 

(first  contacts) ) ) ) 


:ontacts) ) ) ) 


(progn 

(setq  mean-var-list  (compute-mean-covariance  (first 

(setq  lat  (obu-data-field-lat  (obu-contact-spatial-data 

(setq  lng  (obu-data-field-lng  (obu-contact-spatial-data 

(setq  state-cov-old  (make-obu-state-field 

:ref-lat  lat 

:ref-lng  lng 

: state      (first  mean-var-list) 

:covariance  (second  mean-var-list) ) ) 
(setq  time-old  (obu-data-field-obs-time 

(obu-contact-spatial-data  (first 
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i)) 


cse-unc 


) 

(setq  count    (min   (ask  ldbm  contacts-from-head-to-tail) 

(-   (list-length  contacts)    II )t 
(dotimes   (i  count  nil) 

(setq  pos    (+    (obu-track-number-of-contacts  trade)    (-  count)    1 


(setq  etc 
(setq  ctc-sd 
(setq  type 
(setq  obs-time 
(setq  lat 
(setq  lng 
(setq  brg 
(setq  maj 
(setq  min 
(setq  cse 
(setq  cse-unc 
(setq  spd 
(setq  spd-unc 


(nth  pos  contacts) ) 
(obu-contact-spatial-data  etc) ) 
(obu-data-f ield-type     ctc-sd) ) 
(obu-data-f ield-obs-tlme  ctc-sd) ) 


ctc-sd) ) 
ctc-sd) ) 
ctc-sd) ) 
ctc-sd) ) 
ctc-sd) ) 
ctc-sd) ) 
ctc-sd) ) 
ctc-sd) ) 
ctc-sd) ) 


(last  contacts) ) ) 


contacts) ) ) ) ) 


(obu-data-f ield-lat 

(obu-data-f ield-lng 

(obu-data-f ield-brg 

(obu-data-f ield-maj 

(obu-data-f ield-min 

(obu-data-f ield-cse 

(obu-data-f ield-cse-unc 

(obu-data-f leld-spd 

(obu-data-f leld-spd-unc 
(setq  delta-time  (-  obs-time  time-old) ) 
(setq  state-cov-old  (change-tangent  lat  lng 

(ioumotion  delta-time  state-cov-old) ) ) 
(setq  state-cov-ctc  (make-state  type  lat  lng  brg  maj  min  cse 

spd  spd-unc) ) 
(setq  state-cov-old  (center-tangent 

(filter  state-cov-old  state-cov-ctc))) 
(setq  time-old  obs-time) 
) 
(setq  track  (update-track-head  track  uc  state-cov-old 

(obu-contact-spatial-data  (first 

(obu-contact-id  (first  (last 

(setq  track  (update-track-tail  track)) 


(return-from  update-track  track) 


(defobfun  (update-track-head  act-mtst)  (track  uc  state-cov  spatial-data  id) 
(let  ((contacts  (append  (obu-track-contacts  track)  (list  uc) ) ) 
) 
(setq  contacts  (sort  contacts  ♦* (lambda  (x  y) 

(<  (obu-data-field-obs-time  (obu-contact-spatial-data  x) ) 
(obu-data-f ield-obs-tlme  (obu-contact-spatial-data 
y)) )))) 

(setf  (obu-track-number-of-contacts  track)  (♦  (obu-track-number-of-contacts 
track)  1)) 

(setf  (obu-track-head-state-covariance  track)  state-cov   ) 
(setf  (obu-track-head-spatial-data     track)  spatial -data) 
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(setf    (obu-track-head-contact-id  track)    id 

(setf    (obu-track-contacts  track)    contacts 

(return-from  update-track-head  track) 


(defobf 
(let 


y))))) 


un  (update-track-tail  act-mtst)  (track) 

( (nctcs  (obu-track-number-of-contacts  track) ) 
(head-sc  (obu-track-head-state-covariance  track) ) 
(head-sd  (obu-track-head-spatial-data  track) ) 
(head-id  (obu-track-head-contact-id  track) ) 
(tail-sc  (obu-track-tall-state-covariance  track) ) 
(tail-sd  (obu-track-tail-spatial-data  track) ) 
(tail-id  (obu-track-tail-contact-ld  track) ) 
(tail-update-tlme  (obu-data-field-obs-time  tail-sd) ) 
(contacts  (obu-track-contacts  track) ) 

type  obs-tlme  lat  lng  brg  raaj  min  cse  cse-unc  spd  spd-unc 
mean-var-list  new-tail  pos-new-tail  new-tail-sd  delta-time 
state-cov-old  state-cov-ctc  state-cov-new 
) 
(setq  contacts  (sort  contacts  #' (lambda  (x  y) 

(<  (obu-data-field-obs-time  (obu-contact-spatial-data  x)) 
(obu-data-field-obs-time  (obu-contact-spatial-data 

) 

(if  (<-  nctcs  1) 

(progn  (setf  (obu-track-tail-state-covaxlance  track)  head-sc) 
(setf  (obu-track-tail-spatial-data  track)  head-sd) 
(setf  (obu-track-tail-contact-ld  track)  head-id) 
(return-from  update-track-tail  track) 


) 


) 


contacts) ) ) ) 
contacts) ) ) ) 


(if  (<»  nctcs  (ask  ldbm  contacts-frcm-head-to-tall) ) 
(if  (»  tail-id  (obu-contact-id  (first  contacts))) 
(return-from  update-track-tail  track) 
(progn 

(setq  mean-var-list  (ccntpute-mean-covariance  (first  contacts))) 
(setq  lat  (obu-data-field-lat  (obu-contact-spatial-data  (first 


contacts) ) ) 


(setq  lng  (obu-data-field-lng  (obu-contact-spatial-data  (first 

(setf  (obu-track-tall-state-covariance  track)  (make-obu-state-field 

:ref-lat  lat 

:ref-lng  lng 

: state      (first  mean-var-list) 

:covarlance  (second  mean-var-list) ) ) 
(setf  (obu-track-tail-spatial-data  track)  (obu-contact-spatial-data 

(first  contacts) ) ) 
(setf  (obu-track-tail-contact-id  track)  (obu-contact-id  (first 

(return-from  update-track-tail  track) 
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(setq  pos-new-tall  (-  (obu-track-number-of-contacts  track) 

(ask  ldbm  contacts-from-head-to-tail) ) ) 
(setq  new-tall  (nth  pos-new-tall  (obu-track-contacts  track) ) ) 
(setq  new-tall-sd  (obu-contact-spatial-data  new-tall)) 
(setq  type     (obu-data-field-type     new-tall-sd)) 
(setq  obs-time  (obu-data-f ield-obs-time  new-tall-sd) ) 
(setq  lat      (obu-data-field-lat      new-tall-sd) ) 
(setq  lng      (obu-data-field-lng      new-tall-sd)) 
(setq  brg      (obu-data-f leld-brg      new-tall-sd) ) 
(setq  ma j      (obu-data-field-maj      new-tall-sd)) 
(setq  mln      (obu-data-field-mln      new-tail-sd) ) 
(setq  cse      (obu-data-field-cse      new-tail-sd)) 
(setq  cse-unc   (obu-data-field-cse-unc  new-tail-sd)) 
(setq  spd      (obu-data-field-spd     new-tail-sd) ) 
(setq  spd-unc   (obu-data-field-spd-unc  new-tail-sd) ) 
(setq  delta-time  (-  obs-time  tail-update-time) ) 
(setq  state-cov-old  (change-tangent  lat  lng  (ioumotlon  delta-time 

tail-sc) ) ) 

(setq  state-cov-ctc  (make-state  type  lat  lng  brg  naj  min  cse  cse-unc  spd 

spd-unc) ) 

(setq  state-cov-new  (center-tangent  (filter  state-cov-old 

state-cov-ctc) ) ) 

(setf  (obu-track-tail-state-covariance  track)  state-cov-new) 
(setf  (obu-track-tail-spatial-data  track)  new-tail-sd  ) 
(setf  (obu-track-tail-contact-id      track)  (obu-contact-id  new-tail) ) 

(prune-contacts-before-tail  track)   ;  Get  rid  of  contacts  before  tail. 

(return-from  update-track-tail  track) 

) 
) 


(defotfun  (prune-contacts-before-tail  act-mtst)  (track) 
(let  ( (nctcs    (obu-track-number-of-contacts  track)) 
(contacts  (obu-track-contacts         track)) 
(max-ctcs  (ask  ldbm  contacts-from-head-to-tail)) 
tail-position 
) 

(setq  tail-position  (-  nctcs  max-ctcs)) 
(if  (<-  tail-position  0) 

(return-from  prune-contacts-before-tail  nil) 
) 
(dotimes  (1  tail-position  nil) 

(setq  contacts  (remove-element  contacts  0)) 
) 

(setf  (obu-track-number-of-contacts  track)  max-ctcs) 
(setf  (obu-track-contacts         track)  contacts) 
) 
) 


(defobfun  (ioumotion  act-mtst)  (time  state-cov) 
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(let    (Uat  (obu-state-field-ref-lat  state-cov)) 

(lng  (obu-state-field-ref-lng  state-cov) ) 

(mean  (obu-state-field-state  state-cov)) 

(var  (obu-state-f ield-covariance  state-cov) ) 


motion-matrices  phi   f 
vartmp  phitrn 


new-mean  new-vax 


(if 


) 


(-  time  0) 

(return-from  ioumotion  state-cov) 


(setq  motion-matrices  (calc-mot-mat  tijne)  ) 

(setq  phi  (first  motion-matrices)) 

(setq  f   (second  motion-matrices)) 

(setq  new-mean  (matrix-multiply  phi  mean)) 

(setq  vartmp   (matrix-multiply  phi  var  )) 

(setq  phitrn   (matrix-transpose  phi)) 

(setq  new-var  (matrix-add  (matrix-multiply  vartinp  phitrn)  f ) ) 

(setq  phi  new-var) 

(setq  phitrn  (matrix-transpose  new-var)) 

(setq  new-var  (matrix-scalar-multiply  0.5  (matrix-add  phi  phitrn))) 

(return-from  ioumotion  (maxe-obu-state-field 

:ref-lat  lat 

:ref-lng  lng 

: state  new-mean 

:covariance  new-var) 
) 


(defobfun  (calc-mot-mat  act-mtst)  (time) 
(let*  ((invalpha  (/  1  alpha)) 

(invalpha2  (*  invalpha  invalpha)) 
(expalphat  (exponential  (*  -1  alpha  time))) 
(exp2alphat  (*  expalphat  expalphat)) 
(f   (make-array  '(4  4)  : initial-element  0) ) 
(phi  (make-array  '(4  4)  : initial-element  0) ) 


) 

(if  (>  time 
(progn 
(setf 

(setf 
(setf 

(setf 
(setf 
(setf 
(setf 
(setf 

) 

(progn 
(setf 

(setf 
(setf 


0) 

(aref 

<-  I* 

(aref 
(aref 
(*  5 
(aref 
(aref 
(aref 
(aref 
(aref 


(aref 
(-  (* 
(aref 
(aref 


0  0) 
(-  1 

1  1) 

0  2) 


(*  invalpha2  (-  time 


expalphat) )  (« 
(aref  f  0  0) ) 
(•  invalpha2  (^ 


.5  (- 


exp2alphat))) ) 


0) 


(*  invalpha2  (- 
expalphat  1) )  (*  . 

1)  (aref  f  0  0)) 

2)  (*  invalpha2  (♦ 


1  invalpha 
exp2alphat) )>)))) 


0.5  I-  expalphat) 


(aref  f  0  2) ) 
(aref  f  0  2)) 
(aref  f  0  2)) 
(*  0.5  invalpha 
(aref  f  2  2) ) 


(-  1  exp2alphat) ) ) 


(-  time)  (•  Invalpha 

5  (-  exp2alphat  1) ) ) ) ) ) ) 

(-  0.5)  expalphat 


71 


ACT-MTST  Wad,  Apr  25,  1990 


(«  -1  .5  exp; 

>alphat) ) ) ) 

(setf 

(aref  f  2  0) 

(aref  f  0  2) ) 

(setf 

(aref  f  1  3) 

(aref  f  0  2) ) 

(setf 

(aref  f  3  1) 

(aref  f  0  2)) 

(setf 

(aref  f  2  2) 

(*  0.5  invalpha 

(-  exp2 alpha t  1)) ) 

(setf 

(aref  f  3  3) 

(aref  f  2  2) ) 

) 
) 

(setq  f  (matrix-scalar-multiply  (*  sigma  sigma)  f ) ) 

(setq  phi  (identity-matrix  4)) 

(setf  (aref  phi  0  2)  (*  invalpha  (-  1  expalpnat))) 

(setf  (aref  phi  1  3)  (aref  phi  0  2)) 

(setf  (aref  phi  2  2)  expalpnat) 

(setf  (aref  phi  3  3)  (aref  phi  2  2)) 

(return-from  calc-mot-mat  (list  phi  f ) ) 


(defobfun  (exponential  act-mtst)  (x) 
(if  (<  x  -350) 

(return-from  exponential  (exp  -350) ) 
) 
(if  (>  x  350) 

(return-from  exponential  (exp  350)) 
(return-from  exponential  (exp  x) ) 
) 
) 


(defobfun  (make-state  act-mtst)  (type  lat  lng  brg  maj  min  cse  cse-unc  spd  spd-unc) 
(let  ((mean  (make-array   4    : initial -element  0)) 
(var   (make-array  '(4  4)  : initial -element  0)) 
(stationary  (/  (*  sigma  sigma)  (*  2  alpha))) 
(smaj  (max  maj  1))  (smin  (max  min  1)) 
) 

(setf  (aref  mean  2)  (*  spd  (cos  (*  dtor  (-  90  cse))))) 
(setf  (aref  mean  3)  (*  spd  (sin  (*  dtor  (-  90  cse))))) 
(setq  var  (make-var  stationary  brg  smaj  smin)) 
(return-front  make-state  (make-obu-state-field 
:ref-lat  lat 
:ref-lng  lng 
: state  mean 
:covar lance  var)) 
) 
) 


(defobfur.  (make-var  act-mtst)  (stationary  brg  smaj  smin) 
(let*  ((var  (make-array  '(4  4)  : initial-element  01) 
(cost  (cos  (*  dtor  brg) ) ) 
(sinb  (sin  (*  dtor  brg))) 
(cosb2  (*  cosb  cosb) ) 
(sinb2  (*  sinb  sinb) ) 
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(smaj2    (*   0.25  smaj  smajU 
(smin2    (*   0.25  smin  smln) ) 


(setf 

(aref 

var 

0 

0) 

(setf 

(aref 

var 

0 

1) 

(setf 

(aref 

var 

1 

0) 

(setf 

(aref 

var 

1 

1) 

(setf 

(aref 

var 

2 

2) 

(setf 

(aref 

var 

3 

3) 

(return-from  make-var  var) 


(♦    (*   smin2  cosb2)    (*   smaJ2   sinb2) ) ) 

(*    (-  smaj2  smin2)    sinb  cost)) 

(aref  var  0  1)) 

(+  (*  smin2  sinb2)  (*  smaj2  cosb2) ) ) 

stationary) 

stationary) 


(defobfun  (filter  act-mtst)  (state-cov-sol  state-cov-ctc) 


(let 


(Hat 
(lng 
(meanl 
(varl 
(mean2 
(var2 
(infmatl 
(infmat2 
(infvecl 
(infvec2 
(tmpmat 
(tmpvec 
mean  var 


(obu-state-field-ref-lat  state-cov-sol) ) 

(obu-state-field-ref-lng  state-cov-sol) ) 

(obu-state-field-state  state-cov-sol) ) 

(obu-state-field-covariance  state-cov-sol) ) 

(obu-state-field-state  state-cov-ctc) ) 

(obu-state-field-covariance  state-cov-ctc) ) 
(make-array  ' (4  4) ) ) 
4))) 


(make-array  ' (4 

(make-array  4) ) 

(make-array  4) ) 

(make-array  '  (4 

(make-array  4) ) 


4))) 


) 

(setq  infmatl  (matrix- invert  varl)) 

(setq  Infvecl  (matrix-multiply  infmatl  meanl)) 

(setq  infmat2  (matrix-invert  var2) ) 

(setq  infvec2  (matrix-multiply  lnfmat2  mean2)) 

(setq  tinpmat  (matrix-add  infmatl  infmat2) ) 

(setq  var  (matrix-invert  tinpmat)  ) 

(setq  tmpvec  (matrix-add  infvecl  lnfvec2) ) 

(setq  mean  (matrix-multiply  var  tmpvec)) 

(return-from  filter  (make-obu-state-field 

:ref-lat  lat 

:ref-lng  lng 

: state  mean 

:covariance  var) 
) 


(defobfun  (change-tangent  act-mtst)  (newlat  newlng  state-cov) 

(let  ( (oldlat  (obu-state-field-ref-lat    state-cov)) 

folding  (obu-state-field-ref-lng    state-cov) ) 

(mean  (obu-state-field-state     state-cov) ) 

(var  (obu-state-field-covariance  state-cov) ) 

(rot  (make-array  '(4  4)  : initial -element  0)) 

(rottrn  (make-array  "(4  4))) 

(newrmp  (make-array  '(4  4))) 
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(new-mean  (make- array   4    )  ) 

(new-var   (make -array  '(44))) 

latlngllst  rngfcrgllst  cenlat  cenlng  deltabrg  cosb  sinb 

rngl  rng2  rng3  brgl  brg2  brg3 
) 
(setq  rngl  (sqrt  (+  (*  (aref  mean  0)  (aref  mean  0)) 

(*  (aref  mean  1)  (aref  mean  1))))) 
(setq  brgl  (arctan  (aref  mean  1)  (aref  mean  0))) 
(setq  latlnglist  (getlatlng  oldlat  oldlng  rngl  brgl  0)) 
(setq  cenlat  (first  latlngllst)) 
(setq  cenlng  (second  latlngllst)) 

(setq  brgl  (second  (getrngbrg  oldlat  oldlng  cenlat  cenlng  0))) 
(setq  brg2  (second  (getrngbrg  cenlat  cenlng  oldlat  oldlng  0))) 
(setq  rngbrglist  (getrngbrg  newlat  newlng  cenlat  cenlng  0)) 
(setq  rng3  (first  rngbrglist)) 
(setq  brg3  (second  rngbrglist)) 

(setq  brg4  (second  (getrngbrg  cenlat  cenlng  newlat  newlng  0))) 
(setq  deltabrg  (mod  (+  (-  brg2  brgl)  (-  brg3  brg4)  360)  360)) 
(setq  cosb  (cos  (*  dtor  deltabrg))) 
(setq  sinb  (sin  (*  dtor  deltabrg))) 

Move  the  mean  position: 
(setf  (aref  new-mean  0)  (*  rng3  (cos  (*  dtor  (-  90  brg3))))) 
(setf  (aref  new-mean  1)  (*  mg3  (sin  (*  dtor  (-  90  brg3))))) 

Rotate  the  velocity: 
(setf  (aref  new-mean  2)  (♦  <*  cosb  (aref  mean  2))  (•  sinb  (aref  mean  3)))) 
(setf  (aref  new-mean  3)  (-  (*  cosb  (aref  mean  3))  (*  sinb  (aref  mean  2)))) 

Rotate  the  covariance: 

(setf  (aref  rot  0  0)  cosb) 

(setf  (aref  rot  1  1)  cosb) 

(setf  (aref  rot  2  2)  cosb) 

(setf  (aref  rot  3  3)  cosb) 

(setf  (aref  rot  0  1)  sinb) 

(setf  (aref  rot  2  3)  sinb) 

(setf  (aref  rot  1  0)  (-  sinb)) 

(setf  (aref  rot  3  2)  (-  sinb)) 

(setq  newtmp  (matrix-multiply  rot  var) ) 

(setq  rottrn  (matrix-transpose  rot) ) 

(setq  new-var  (matrix-multiply  newtmp  rottrn) ) 

(setq  rot  new-var) 

(setq  rottrn  (matrix-transpose  new-var) ) 

(setq  new-var  (matrix-scalar-multiply  0.5  (matrix-add  rot  rottrn))) 

(return-from  change-tangent  (make-obu-state-field 

: ref-lat    newlat 

:ref-lng    newlng 

: state     new-mean 

: covariance  new-var) 
) 
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(defobfun   (center-tangent  act-mtst)    (state-cov) 

(let    (Uat      (obu-state-field-ref-lat         state-cov)) 

(lng     (obu-state-field-ref-lng        state-cov)) 

(mean  <obu-state-fleld-state     state-cov) ) 

(vaz  (obu-state-field-covaxlance  atate-eov) ) 

(rot   (make-array  '(4  4)  : Initial -element  0)) 

(rottrn  (make-array  '(4  4))) 

(newtmp  (make-array  '(4  4))) 

(new-mean  (make-array   4   ) ) 

(new-var   (make-array  '(4  4))) 

latlnglist  rngbrglist  newlat  newlng  newrng  newbrg 

rngl  brgl  cost  sinb 
) 
(setq  rngl  (sqrt  (♦  (*  (aref  mean  0)  (aref  mean  0)) 

(*  (aref  mean  1)  (aref  mean  1))))) 
(if  (<  rngl  0.01) 

(return-froni  center -tangent  state-cov) 
) 

(setq  brgl  (arctan  (aref  mean  1)  (aref  mean  0))) 
(setq  latlnglist  (getlatlng  lat  lng  rngl  brgl  0)) 
(setq  newlat  (first  latlnglist)) 
(setq  newlng  (second  latlnglist)) 

(setq  rngbrglist  (getrngbrg  newlat  newlng  lat  lng  0)) 
(setq  newrng  (first  rngbrglist) ) 
(setq  newbrg  (mod  (♦  (second  rngbrglist)  180)  360)) 
(setq  deltabrg  (-  newbrg  brgl)) 
(setq  cosb  (cos  (*  dtor  deltabrg) ) ) 
(setq  sinb  (sin  (*  dtor  deltabrg) ) ) 

Move  the  mean  position: 
(setf  (aref  new-mean  0)  0) 
(setf  (aref  new-mean  1)  0) 


Rotate  the  velocity: 


(setf 

(aref  new-mean  2)  <♦  (*  cosb  (aref  mean  2) )  ( 

(setf 

(aref  new-mean  3)  (-  (*  cosb  (aref  mean  3) )  ( 

Rotate  the  cover lance: 

(setf 

(aref  rot  0  0)  cosb) 

(setf 

(aref  rot  1  1)  cosb) 

(setf 

(aref  rot  2  2)  cosb) 

(setf 

(aref  rot  3  3)  cosb) 

(setf 

(aref  rot  0  1)  sinb) 

(setf 

(aref  rot  2  3)  sinb) 

(setf 

(aref  rot  1  0)  (-  sinb)) 

sinb  (aref  mean  3) ) ) ) 
sinb  (aref  mean  2)))) 


(setf  (aref  rot  3  2)  (-  sinb)) 

(setq  newtmp  (matrix-multiply  rot  vax) ) 

(setq  rottrn  (matrix-transpose  rot) ) 

(setq  new-var  (matrix-multiply  newtmp  rottrn) ) 

(setq  rot  new-var) 

(setq  rottrn   (matrix-transpose  new-var)) 

(setq  new-var   (matrix-scalar-ttultiply  C.5   (matrix-add  rot  rottrn))) 

(return-from  center-tangent    (make-obu-state-field 

:ref-lat         newlat 

:  ref  -lng         newlng 
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:  state  nwMiii 

:covarlance  nev-var) 
) 

) 


(defobfun    (arctan  act-mtst)    (y  x) 
(let    (a 
) 
(If    (and   (-  y  0)    (•  x  0) ) 

(return- from  arctan  0) 
) 
(if   (-  y  0) 

(progn   (if    (<  x  0) 

(retum-from  arctan  270) 
( return- f ran  arctan     90) 
) 
) 
) 
(if   (<  y  0) 

(setq  a  (♦  (/  (at an  (/  x  y) )  dtor)  180)) 
(setq  a    (/  (atan  (/  x  y) )  dtor)     ) 
) 
(If  (>  a  0) 

(return-frcm  arctan  a) 
(return-froa  arctan  <♦  a  360) ) 
) 
) 
) 


This  procedure  converts  a  cover  lance  matrix  (cov)  into  an 
ellipse.  The  ellipse  is  returned  In  the  form: 

(  Orientation  .  Semi-major  axis  ,  Semi -minor  axis  ) . 


(defobfun  (cov-to-elllpse  act-mtst)  (cov) 
(let  ((a  (aref  cov  0  0)) 
(b  (aref  cov  0  1) ) 
(c  (aref  cov  ID) 

brg  cosb  slnb  cosb2  slnb2  slncosZb  smjr  smnr 
) 
(If  (not  (>  (-  (•  a  o  (•  b  b))  0)) 

(return-f rom  cov-to-elllpse  (list  0  (expt  10  9)  (expt  10  9) ) ) 
) 
(if  (-  b  0) 

(If  (not  (>  a  c>> 

(return-frcm  cov-to-elllpse  (list  0 

(return-from  cov-to-elllpse  (list  90 


(•  2 

(sqrt 

(abs  0 ) ) 

C  2 

(sqrt 

(abs  a) ) ) ) ) 

C  2 

(sqrt 

(abs  a))) 

C  2 

(sqrt 

(abs  odd 

) 

(If 

(>  b 
(If 

) 

0) 
(-  a  c) 

( return - 

■from 

cov 

■to-elllpse 

(list 

45 

C 

2 
2 

) 

(if 

«  b 

0) 

(sqrt  (♦  a  b)>) 
(sqrt  (-  a  b) ) ) ) ) 


(If  (-  a  c) 

(return-from.  cov-to-elllpse  (list  135  (*  2  (sqrt  (-  a  b))> 

(•  2  (sqrt  (♦  a  b))))) 

) 
) 

(setq  brq   (/    (arctan   (-  c  a)    (*   2  bl)    2)1 
(setq  slnb    (sin    (*  dtor  brq))) 
(setq  cosb    (cos    (*  dtor  brq))) 
(setq  slnb2    (*  slnb  slnb)) 
(setq  cosb 2    (*  cosb  cosb)) 
(setq  slncos2b    (•  2  b  slnb  cosb)) 

(setq  smjr    (•  2    (sqrt    (♦    (•  a  slnb2)    slncos2b   (•  c  cosb2))))) 
(setq  smnr    (•  2    (sqrt    (♦    (•  a  cosb2)    (•  -1  slncos2b)    (*  c  sinb2)  ID) 
(return-from  cov-to-elllpse    (list  brq  smjr  smnr)) 
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.-lisp  matrix  manipulation  routines 
;4/28/89--s.c.  lent 

;  these  routines  operate  on  three (3)  'standard  type'  lisp  objects  which  are 
;  subsets  of  lisp  object  types: 

1.  scalar s,  i.e.,  lisp  type  'number' 

2.  vectors.  I.e..  lisp  type  'simple-vector'  with 

length  >  1 

3.  matrices,  i.e.,  lisp  type  'array'  with 
;  rank  -  2  and 

each  dimension  >  1 

;  all  routines  accept  input,  and  produce  output  of  these  types 

;  unacceptable  inputs  produce  'MIL'  as  output,   there  are  two (2)  cases: 

1.  'non-standard'  inputs 

2.  'standard'  but  mismatched  inputs,  e.g., 
the  Inverse  of  matrix (m  x  n),  where  m  o  n. 

.•routines  consist  of: 
;  1 .  mx-mult 

;  2 .  mx-add 

r  3.  sc-mx-mult 

4.  mx-trans 

5.  mx-ident 

6.  mx-inverse 

these  are  built  using  the  following  'sub-routines': 

1.  matrlxp 

2.  array-row-el 

3.  array-col -el 

4 .  array-row-dim 

5.  array-col-dlm 

6.  can-mx-mult 

7 .  can-mx-add 

8.  sqr-mx 

9.  mx-f actor 
10.  mx-subst 


(defun  matrix-multiply  (a  b) 
(if  (not  (can-mx-mult  a  b) ) 

(return-from  matrix-multiply  NIL) 
) 

(let  ((c)  (ij-list)) 
(setq  c  (make-array  (remove  1  (list  (array-row-dim  a) (array-col-dim  b) ) ) 
) 
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) 

(do  ((indexl  0  (+  indexl  1))) 

((>  indexl  (-  (array-col-dim  b)  1))) 

(if  (>  (array-col-dim  b)  1)  (setq  ij-llst  (list  indexl)) 

(setq  ij-list    ()) 
) 
(do    <<iridex2  0    (+  index2  1))) 

( (>  index2    (-    (array-row-dim  a)    1))) 

(if  (>  (array-row-dim  a)  1)  (setq  ij-list  (cons  index2  ij-list)) 

) 

(setf  (apply  #'aref  c  ij-list)  0) 

(do    <(index3  0    (+  index3  1))) 

((>  index3  (-  (array-dimension  b  0)  1))) 

(setf  (apply  i'aref  c  ij-list)  (+  (apply  i'aref  c  ij-list) 

(*  (axray-col-el  a  (list  index2 


index3) ) 
indexl)) 


) 


(array-row-el  b  (list  index3 


) 


) 

(if  (>  (array-row-dim  a)  1)  (setq  ij-llst  (cdr  ij-list)) 

) 


) 

(if  (equal  (array-dimensions  c)  0  ) 
(aref  c) 

c 
) 


(defun  matrix-scalar-multiply  (sc  mx) 
(if  (not  (and  (numberp  sc) (matrixp  mx))) 

(return-from  matrix-scalar-multiply  KID 
) 

(let  ( (ij-list) (mx-new  (make-array  (array-dimensions  mx) )) ) 
(do  ((indexl  0  (+  indexl  1))) 

( (>  indexl  (-  (array-col-dlm  mx)  1) ) ) 

(if  (>  (array-col-dim  mx)  1)  (setq  ij-list  (list  indexl)) 

(setq  ij-list  ()) 


) 
(do 


((lndex2  0  <+  index2  1))) 

((>  lndex2  (-  (array-dimension  mx  0)  1))) 

(if  (>  (array-dimension  mx  0)  1)  (setq  ij-list  (cons 


) 

(setf  (apply  I'aref  mx-new  ij-list) 

(setq  ij-list  (cdr  ij-list)) 


index2  ij-list) ) 
(apply  #'aref  mx  ij-list)  so) 
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(defun  matrix-transpose  (mx) 
(if  (not  (matrixp  mx) ) 

(return-from  matrix-transpose  NIL) 
) 

(let  ((ij-list)  (mx-new  (make-array  (reverse  (axray-dimenslons  mx)  )  )  )  i 
(do  ((indexl  0  (+  indexl  1))) 

( (>  indexl  (-  (array-col -dim  mx)  1))) 

(if  (>  ( array-col -dim  mx)  1)  (setq  lj-llst  (list  indexl)) 

(setq  lj-llst  ()) 
) 
(do  <(lndex2  0  (♦  lndex2  1))) 

((>  index2  (-  (array-dimension  mx  0)  1))) 

(if  (>  (array-dimension  mx  0)  1)  (setq  lj-list  (cons  lndex2  ij-llst)) 
) 

(setf  (apply  t'aref  mx-new  (reverse  ij-list))  (apply  #'aref  mx  ij-list)) 
(setq  ij-list  (cdr  ij-list)) 
) 
) 

mx-new 
) 
) 


(defun  matrix-add  (mxl  mx2) 
(if  (not  (can-mx-add  mxl  mx2)) 

(return-from  matrix-add  NIL) 
) 

(let  ( (lj-llst) (mx-new  (maXe-array  (array-dimensions  mxl) )) ) 
(do  ((indexl  0  (+  Indexl  ID) 

((>  indexl  (-  (array-col -dim  mxl)  1))) 

(if  (>  (array-col -dim  mxl)  1)  (setq  ij-list  (list  indexl)) 

(setq  ij-llst  ()) 
) 
(do  ((index:  0  (♦  index2  1))) 

((>  lndex2  (-  (array-dimension  mxl  0)  1))) 

(if  (>  (array-dimension  mxl  0)  1)  (setq  ij-list  (cons  index2  ij-list)) 

) 

(setf  (apply  i'aref  mx-new  ij-list)  (♦  (apply  I'aref  mxl  ij-llst) 

(apply  i'aref  mx2  ij-list))) 
(setq  ij-list  (cdr  ij-llst)) 
) 
) 

mx-new 
) 
) 


(defun  identity-matrix  (m) 
(if  (not  (and  (lntegerp  m) (>  m  1))) 

(return-from  identity-matrix  NIL) 
) 

(let  ( (ij-list)  (mx-new  (maXe-array  (list  ma)))) 
(do  ((indexl  0  <♦  indexl  1))) 
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((>  lndexl  (-ml))) 
(setq  lj-llst  (list  lndexl)) 
(do  <<lndex2  0  (+  lndex2  1))) 
((>  index2  (-  m  1))) 
(setq  lj-llst  (cons  lndex2  ij-llst)) 

(setf  (apply  #'ar*f  u-im  lj-llst)  (If  (-  lndexl  lndex2)  1  0)) 
(setq  lj-llst  (cdr  lj-llst)) 
) 
) 

mx-new 
) 
) 


(defan  matrix-Invert  (a) 
(if  (not  (sqr-mx  a) ) 

(return-from  matrix -invert  MIL) 
) 

(let*  ( (n  (array-dimension  a  0)) 
(temp    (mx-f actor  a  n) ) 
(w      (car  temp) ) 
(ipivot  (cadr  temp) ) 
) 
(if  (-  0  (caddr  temp) ) 

(return-from  matrix-invert  NIL) 
) 

(mx-subst  w  n  ipivot) 
) 
) 


(defun  matrixp  (a) 
(if  (arrayp  a) 

(cond  ( (simple-vector-p  a  ) 

(if  (>  (array-dimension  a  0)  0) 
T 

NIL 
) 
) 

( (equal (array-rank  a)  2) 
(if  (and  (>  (array-dimension  a  0)  0) 
(>  (array -dimension  a  1)  0) 
) 

T 

NIL 
) 
) 

(T  NIL) 
) 
NIL 
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(defun  array-row-el   (a  1) 

(apply  *'aref  a  (if  (-  (array-rank  a)  1) 

(Use  (car  i)  ) 
1 
) 
) 
) 


(defun  array-col -el   (a  i) 

(apply  I'aref  a  (if  (•  (array-rank  a)  1) 
(cdr  1) 

1 

) 
) 
) 


(defun  array-row-dim  (a) 
(If  (-  (array-rank  a)  1) 
1 

(array-dimension  a  0) 
) 
) 


(defun  array-col -dim  (a) 
(if  (-  (array-rank  a)  1) 
1 

(array-dimension  a  (-  (array-rank  a)  1)) 
) 
) 


(defun  can-mx-mult  (a  b> 

(if  (and  (macrixp  a) (matrixp  b) ) 

(if  (-  (array-dimension  a  (-  (array-rank  a)  1) ) 
(array-dimension  b  0) 
) 

T 

NIL 
) 
NIL 


(defun  ean-tfuc-add  (a  b) 

(if  (and  (matrixp  a) (matrixp  b) ) 

(if  (and  (-  (array-dimension  a  0) 
(array-dimension  b  0) 
) 


81 


Matrix  (unctions  Wad,  Hay  3.  ItM 


(-  (array-dimension  a  (-  (array-rank  a)  1)) 
(array-dimension  b  (-  (array-rank  b)  1)) 
) 
) 

T 

NIL 
) 
NIL 

) 


(defun  sqr-mx    (a) 
(If    (matnxp  a) 

(if    (-   (array-rank  a)    2) 

(If  (-  (array-dimension  a  0) (array-dimension  a  1) ) 

T 

NIL 

) 
) 
NIL 


(defun  mx-factor  (mx  n) 
(let  ( (awlkod) (colmax i (ratio) (rowmax) (temp)  (lstar) (iflag  1) 
(iplvot  (make-array  n) ) 
(d      (make-array  n) ) 

(w      (make-array  (array-dimensions  mx) ) ) 
) 
(if  (<-  n  0) 

(return-from  mx-f actor  nil) 
) 


(do  ((10  (+11)) 
(rowmax  0  0) 
) 

((-in)) 

(setf  (aref  Iplvot  1)  1) 
(do  ((j  0  (♦  j  1))) 
((-  j  nil 

(setf  (aref  w  l  j)  (aref  mx  1  j)) 
(setq  rowmax  (max  rowmax  (abs  (aref  w  1  j)))) 
I 
(when  (-  rowmax  0) 

(setq  lflag  0) 
(setq  rowmax  1) 
) 
(setf  (aref  d  1)  rowmax) 


(do  <<k  0  (♦  k  1))) 
((-  k  (-  n  1))) 
(setq  colmax  (abs  (/  (aref  w  k  k) (aref  d  k) ) ) ) 
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3)))) 


(setq  lstar  k) 
(do  ((  l  (♦  k  1)  (♦  i  1))) 
((•  1  n)) 

(setq  awlkod  (abs  (/  (aref  w  1  k) (aref  d  1)))) 
(when  (>  awlkod  colaax) 

(setq  colmax   awlkod) 
(setq  lstar  1) 
) 
) 

(If  (-  colmax  0) 
(setq  lflag  0) 
(progn 

(If  (>  lstar  k) 
(progn 

(setq  lflag    (*   -1   If  lag) I 
(setq  1    (aref  lplvot   lstar) ) 
(setf    (aref  lplvot  lstar)    (aref  lplvot  k) ) 
(setf    (aref  lplvot  k)    1) 
(setq  temp   (aref  d  lstar) ) 
(setf    (aref  d  lstar)    (aref  d  k) ) 
(setf    (aref  d  k)    temp) 
(do  <(}  0   (♦  3  1))) 
((-  J  ni) 

(setq  temp    (aref  w  lstar   }) ) 
(setf    (aref  w  lstar  j) (aref  w  k  J)) 
(setf    (aref  «  k  J|    temp) 
) 
) 
) 

(do    ((1    ckll    IM1III 
((-  1  n)) 

(setf  (aref  w  1  k)  </  (aref  w  1  k)  (aref  w  k  k) ) ) 
(setq  ratio  (aref  w  1  k) ) 
(do  ((  3  (♦  k  1)  (♦  Jl))) 
((-  3  nl) 
(setf  (aref  »l  ])  (-  (aref  •  1  )i  (*  ratio  (aref  w  k 


) 


) 


) 


) 

(If  (-  (aref  w  (-  n  1)<-  n  D)  0) 

(setq  lflag  0) 
) 
(list  w  lplvot  lflag) 


(defun  mx-subst  (w  n  lplvot) 
(let  ( 

(b      (make-array  n  : Initial -element  0)) 
(ainv   (make-array  (list  n  n) ) ) 
dp) 
) 

(do  ((col  0  (♦  col  1))) 
((-  col  n)) 
(setf  (aref  b  col)  1) 

(setq  ip  (aref  lplvot  0) ) 

(setf  (aref  ainv  0  col) (aref  b  ip) ) 

(do  111  1  (*  1  1)1 
(sum  0  0) 
) 

((-  1  n)) 
(do  ((3  0  (♦  3  !>>> 

(setq  sun  (♦  sub  (•  (aref  w  1  3>  <*"*  *lnv  1   »1)))J 
) 

(setq  lp  (aref  lplvot  1>> 
(setf  (aref  ainv  i  col  )  (-  (aref  b  lp)  sum)) 

'(setf  (aref  ainv  <-  n  1)  col)  (/  (aref  ainv  (-  n  1)  col)  (aref  w  (-  n 


.)  (-  n  1) ) ) ) 


(do  ((  1  (-  n  2)  (-11)) 
(  sum  0  0) 
) 

((-  1  -1)1 

(do  ((  j  (♦  i  1)  (♦  3  I"' 
((-  3  n)) 
(setq  sum  (♦  sum  (*  (aref  ainv  3  col) (aref  w  1  3)))) 

(setf  (aref  ainv  1  col)  (/  (-  (aref  ainv  1  col)  sum) (aref  w  i  ill: 
) 


) 
ainv 


(setf  (aref  b  col)  0) 
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TRACK  Utility  Procedures: 


(defun  getLatLng  (latl  lngl  dlst  brq  flag) 
: flag  ■  0  great  circle  calc 
-  1  rhumbllne 
brg  in  degrees. 

(let*  (Uat2  latl)  Ung2  lngl)  (mg  (rem  </  dlst  60.1077)  360)) 
(cbrg  (cosLisp  brg))  (latlt  latl)  (clatlt)  (slatl) 
(srng  nil)  (crng  nil)  (ding  nil)  (temp  nil)  (clat2  nil) 
) 
(if  (-  latl  90)  (setq  latlt  89.9)) 
(if  (-  latl  -90)  (setq  latlt  -89.9)) 
(setq  clatl  (cosLisp  latlt)) 
(if  (-  flag  0)      ;  great  circle  calculation 
(progn  (if  (>  rng  180) 

(progn  (setq  rng  (-  360  rng)) 

(aetq  brg  (rea  (♦  180  brg)  360) ) 
(setq  cbrg  (cosLisp  brg) > 
) 
) 
(if  (-  clatl  0) 

(progn  (if  (-  (sinLisp  brg)  1) 

(Mtq  lat2  (-  latl  rng)) 
(setq  lat2  (♦  latl  rng)) 
) 

(setq  lng2  lngl) 
) 

(progn  (aetq  slatl  (ainLisp  latlt)) 
(setq  srng  (ainLisp  rng) ) 
(aetq  crng  (cosLisp  rng) ) 
(aetq  lat2  (-  90 

(acosLlap  (fnAbs  I 


(setq  clat2  (cosLisp  lat2)) 
(setq  ding  (acoaLisp  (fnAbs 


clatl  srng  cbrg) 
slatl  crng) ) > ) ) ) 


(/ 


(-  crng  (•  slatl 

(sinLisp  lat2) ) i 
(•  clatl  clat2) ) i ) ) 


(If  «  (sinLisp  brg)  0) 

(aetq  lng2  (fnCor  (-  lngl  ding) ) ) 
(setq  lng2  (fnCor  <♦  lngl  ding) ) ) 


) 


) 


) 


) 


180 


;lf  clatl-0 
;  progn  rng  > 
)  .-if  flag  -  gc 

(if  (-  flag  1)   ;rhupjblina  calc 
(progn  (if  (or  (•  brg  90) 
(-  brg  270) ) 
(progn  (setq  lat2  latlt) 
(if  (-  brg  270) 

(aetq  lng2  (fnCor  (-  lngl 

(/  rng  clatl) ) ) ) 
(setq  lng2  (fnCor  (♦  lngl 


) 


(/  rng  clatl)))) 


)    .-end  then  construct  brg*>90  or  brg-270 
(progn   (setq  lat2    (♦   latlt    (*  rng  cbrg) ) ) 
(if   (»  (aba  lat2)   90) 
(progn   (if   (>  lat2  0) 

(setq  lat2  90) 
(setq  lat2  -90) 


) 

(aetq  lng2  0) 
) 
(progn  (setq  temp  (log 


(/ 


(tanLisp 
(tanLisp 
(setq  lng2  (fnCor  (♦  lngl  (• 


(-  45 
(-  «5 


latlt) ) ) 
lat2) ) ) ) )  ) 


) 


GLB . radToDeg  temp 
(tanLisp  brgi ) ) ) ) 


) 


)    .-end  else  construct  brg-90  or  brg»270 


) 


) 


)  ;  if  flag-1 
(list  lat2  lng2) 
)   ;  let 
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isefun  getRngSrg  llacl  lngl  lac2  lng2  flag) 
;  flag  -  0  great  clrcla 

-  1  rhunbllne 
(lac*  <<rng  0)  (brg  0)  (dlac  (aba<-  lacl  lat2))) 

(ding  (abs(fnCor(-  lngl  lng2))))  (x)  (y)  (latlt)  (clat2  (cosLisp  lat2  )) 
) 
(If  (»  (♦  dlat  ding)  0.0001) 
(progn  (If  (-  flag  0) 

(progn  (if  (»  ding  01 

(progn  (if  (<  latl  lat2) 
(satq  brg  0) 
isatq  brg  180) 
) 

(satq  mg  dlat) 
) 

(progn  (if  (•  tabs  lacl)  90) 
(secq  laclC  (•  89.9 

(/  lacl  (abs  lacl) ) ) i 
(•acq  lacic  lacl) 
) 

(sacq  claci  (cosLisp  laclC)) 
(sacq  slacl  (SlnLisp  latl) ) 
(sacq  slac2  (slnLisp  lac2)) 
(satq  rng  (acosLisp  (fnA£>s<  +  (•  slatl  slat2) 

(*  datl  clat2 

(satq  brg  (acosLlsp  (fnAbst/  (-  slat2 

(•  slatl  (cosLisp 


(cosLisp  ding) ) ) ) ) ) 


mg))) 
rng) ) ) ) ) ) 

datalina  crossing 

;  is  this  progn  necessary 


(If  (or  (>-  lng2  -90) 
«-  lngl  90)) 


(•  clacl  (slnLisp 


;othar  than  east  or  wast 


(progn  (If  (<  lng2  lngl) 

(satq  brg  (-  360  brg) ) 

(progn  (if  (and  (>  lng2  90) 

(<  lngl  -90)) 
(satq  brg  (-  360  brg) ) ) 


)) 


)) 


)> 
))  ;if  flag-0 
(if  (-  flag  1) 
(progn  (If  (or  <-  (abs  latl)  90) 
(-  (abs  lat2)  90)) 
(progn  (sacq  rng  dlat) 

(if  (>  lac2  lacl) 
(satq  brg  0) 
(sacq  brg  180)) 


r North  or  South  pola 


(progn  (««tq  x  (log  (/  (tanUsp  (-  45  (•  0.5  lacl))) 
iprcvn     h  ~*  (tanU>J  (.  45  {.  0.5  UC2)))))) 

(sacq  x  (•  x  GLB . radToOag) ) 
(satq  y  (fnCor  (-  lng2  lngl))) 
(satq  brg  (fnArctan  x  y)) 
(if  (or  (-  brg  90) 

(-  brg  270)) 

(satq  rng  (•  clac2  ding) ) 

(sacq  rng  (abs  (/  dlat  (cosLisp  brg) ) ) ) 


)) 


)) 


(secq  rng  (•  60.1077  rng)) 
(Use  rng  brg) 
)  .lac 
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Detection  fUport. r  Wed,  Mir  21 , 1  WO 


(setq  detect ionReporter    (kindof  generalPiatform) l 


POS+VTI.      ) 

POS-pos  only. 

POS*VEL-pos   *  vel 

10                 ) 

ret» 

3                 ) 

lets 

30                ) 

degrees 

(defobfun    (exist  detect lonReporter)    (lnlt-llst) 
(usual-exist   lmt-Ilst) 

(have  'report -type 

(have  ' positional -unc 

(have  'speed-unc 

(have  'course-unc 

(have   'expRepDelay     5) 
(have    'FAR  0.9) 

(have   'startTrackp     1) 

(setq     dataNamelist    (append  dataNajneLlst 
1  (:  reports 

(report-type  positional -unc  speed-unc  course-unc 
expRepDelay  FAR  startTrackp) ) ) ) 
(setq     dataValuelist    (append  dataValuelist 
(:  reports 

(.report-type  .positional -unc  , speed-unc  ,  course-unc 
, expRepDelay  , FAR  , startTrackp) i ) ) 
(setq     data?ypeJ.ist    (append  dataTypelist 

' (: reports    (data  data  data  data  data  data  data) ) ) ) 
(setq     dataTextlist    (append  dataTextList 

•CreportS    ("REPCRT  TYPE"   "POSITION  ERROR"    "SPEED  ERROR" 
"COURSE  ERROR*  "Exp  Processing  Delay" 
"False  Aianr.  Rate"   "Track  Initiaticr.")  )  )  ) 
(setq     dataTeirplateText    (append  dataTentplateText 
' ( : reports 

(l-POS  or  POS-VEI."   "[POS   ,    POS-VEI.!") 
("Nautical  Miles"  "(0  -  1000!") 
("Knots"  "[0  -  30i") 

("Degrees"  "(0  -  180]") 

("Hours"  "(0  -  100]") 

("FA  per  Day"  "(0  -  unlimited! ") 

("Single  Contact  to  Track"  "0  -  no,    1  - 

yes") ) ) ) ) 

(let*    ( (rt-designators     • (report s-rb  )) 

(rs-values  (list  durrr.yRBdata) ) 

(rs-names  '  (duwr.yRBdata) ) 

) 

(update-rb-lists  rb-designators  rb-values  rb-namesl 
) 
) 


(defocfun    (setForStart  detect ionReporter)    0 
(usual-setForStart) 
(have    '  cotr.p-jtedTAR  FAR) 
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Detection  Reporter 


•r  21, II 


(defobfun    (getReportDelay  detectionReporter)     (toptional   taxgetType) 
(return-from  getReportDelay 

( exponent lalDraw  expRepOelay) 


) 


) 


(defobfun    (make-obu-report  detect ionfleporter)     (time  target-parma  soptionai   HFDFp) 
(let*    ((cse  (platformStateParms-cse  target -panr.s )  ) 

(spd  (platfcmStateParrrs-spd  taxget-parr.sl) 

(la*  (platforxStateParms-lat  taxget-parms I ) 

(lng  (platformStateParms-lng  target -pains ) ) 

(hgt  <platformStatePanns-hgt  target-perms) ) 

(target Name    (platf  or-tStateParns-platfonnNaae  target-parns) ) 
(cat  (maxe-oDu- category   :pinnedp  nil 

:lockedp  nil 
: cluster-count  0 
: as soc- count        0) ) 
(sensor-field   (make-obu-sensor 

: name    (prlnc-to-strlng  name) 

: single-report-to-trackp    (not    (equalp  startTraekp  0)))) 
spatial -data-field 
) 


3€0) 
)) 


(if  (equalp  report-type  'POS-VEL) 

(setq  cse  (mod   <♦  cse  (-  (random  (•  2.0  course-unc) )  course-unc) ) 

spd  (max  0  (♦  spd  (-  (random  (•  2.0  speed-unc))  speed-unc) ) 


) 


(setq  cse  nil 
spd  nil) 


(multlple-value-setq  (lat  lng) (contectPosltionDraw  iat  lng 
positional -unci ) 

(setq  spatial-data-field     (make-obu-data-field 


:type 

report -type 

:obs-tiae 

time 

:lat 

lat 

ring 

lng 

:brg 

0 

:maj 

positional -unc 

:mln 

posltional-unc 

:cse 

cse 

: cse-unc 

course-unc 

:spd 

spd 

: spd-unc 

speed-unc 

(return- from  make-obu-report  (make-obu-contact 

:id 


nil 
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O»t»etlon  Reporter  Wed.  Mir  21 ,  1M0 


: receipt -tine  time 
:track-eaaoclation  targetNane 

: aenaor  eenaor-f ield 

: categor i zat  ion  cat 

:apatlal-data  apatlal -data-field 

.-altitude  hgt 

:BFDFp  HFDFp 
) 


(defobfun  (maXe-false-obu-report  detectionReporter) (lat  lng  uncertainty  ctime 
targetType ) 

(let  ( (f alaeTargetState  (getFalaetargetState  lat  lng  uncertainty  targetType ) ) ) 
(maXe-obu-report  ctime  faiaeTargetState) 

) 
) 


(defobfun  (getFalaetargetState  detectionReporter)  (lat  lng  uncertainty  targetType) 
(mul t  ipl e - val ue-b  i  nd 
(lat  lng) 

(contactPositionCraw  lat  lng  uncertainty) 
(returr.-from  getFalaeTargetState 
(maJce-platfonnStatePazma 
:lat  lat 
:lng  lng 
:hgt  0 

:cae  (•  (random  1.0)  360) 
:spd  10 

:platfermNaine  'FA 
: targetType  targetType 
) 
) 
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Structures  Thu.  %  10, 1M0 


(def struct  coast -event 

tune 

object 

procedure 

data 

updateList 
) 


(defstruct  coast-message 

send-time 

receipt-time 

type 

sender 

content 

size  :  MBs 

transmission-path 

transmission-count 
) 


Time  message  sent. 

Time  message  received  at  last  node. 

E.g.  'Detect-msg 


(defstruct  obu-contact 

id 

receipt-time 

track-association 

sensor 

categorization 

spatial -data 

altitude 

HFDFp 
) 

(defstruct  obu-track 

Id 

number-of -contact s 

head- state-covar lance 

head-contact -Id 

head-spatial -data 

tai  1  -  state-  cowl  ance 

tail -contact -id 

tall-spatial -data 

altitude 

(contacts  nil  :type  list) 
) 

(defstruct  obu-data-field 
type 

obs-time 
lat 
lng 
brg 
ma  J 

mln 

cse 

cse-unc 
spd 
spd-unc 

) 

(defstruct  obu-state-field 

ref-lat 

ref-lng 

(state     nil  -.type  (array  float    4  )) 

(covariance  nil  :type  (array  float  '(4  4))) 
) 

(defstruct  obu-sensor 
name 
single-report -to-trackp 

) 

(defstruct  obu-category 

pinnedp 

lockedp 

assoc-count 

cluster-count 

ambiguity-list 

ambiguity-resolution 

ambiguity -resolution-time 
) 
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Probability  function*  Frt,  Mar  23, 1990 


(defun  contactPositionDraw  (lat  lng  rng) 

;  a  draw  {rain  a  circular  normal  —  rng  Is  Interpreted  as  a  2  sigma  confidence 
;  Interval 

(let  ( (randomRange  (mln  {*  (/  rng  2)  (sqrt  (*  -2  (log  (random  1.0))))) 
(•  4  rng) ) ) 
(random/ingle  (randonCse) ) 
) 

(return-from  contactPositionDraw 
(values-list 

(getLatLng  lat  lng  randomRange  randomAngle 
•greatCircle*) ) ) 
) 

) 


SPA  FUNCTIONS 


(defun  expandAOUpos  (vel  tau  timeLate) 
(let  ((ratio  (/  timeLate  tau)) 

(factor  (•  2  pi  vel  vel  tau  tau))) 
(*  factor  (+  ratio  -1  (exp  <•  -1  ratio) ) ) ) 
) 
) 


(defun  expandAOUpos vel  (vel  tau  timeLate) 

(let*  ((ratio  (/  timeLate  tau)) 

(factorl  (*  2  pi  vel  vel  tau  tau)) 

(factor2  (*  2   (exp  (*  -1  ratio)))) 

(factor3  (•  C.5  (exp  (•  -2  ratio)))) 
) 

(*  factorl  (♦  ratio  factor2  -1.5  (*  -1  factor3) ) ) 


) 
) 


(defun  getRandomPointlnRegion  (region) 

(let*  ( (regionParms  (getRegionParms  region) ) 

(minLat  (region-parms-mlnLat  regionParms) ) 
(maxLat  (region-parms-maxLat  regionParms) ) 
(minLng  (reglon-parms-minLng  regionParms) ) 
(maxLng  (region-parms-maxLng  regionParms)) 
(ckStuf f  (reglon-parms-checkpts  regionParms) ) 
(thePoint  nil) 
) 

(if  (>  mining  maxlng)  (setq  mining  (-  mining  360))) 
(do  ((insldep  nil)) 
(insldep) 
(setq  thePoint  (list 

(♦  mlnlat  (*  (-  maxlat  minlat)  (random  1.0))) 

(+  mining  (*  (-  maxlng  mining)  (random  1.0))))) 
(setq  insldep  (checlclf  Inside  thePoint  ckStuff ) ) 
) 
(return-from  getRandomPointlnRegion  thePoint) 

) 
) 
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APPENDIX  C:    ITERATIONS  TO  ACHIEVE  STEADY  STATE  GAIN  VALUES 


The  model  in  this  appendix  calculates  the  percent  differece  between  succesive  iterated  values  of  Kalman  gain  K 
These  values  are  tabulated  for  various  ranges  of  interarrival  time  and  sensor  AOU  size. 


0     =    0..25 


A     =    50 


tOO        a   -  .25 


=  fib 


2-or 


rl    =    25  V 


vlQ    -  rl  v2Q   ^0        v30    -  rt  <t>2   =  -•[  I  - exp[-crA]]        *3   -exp[-a-A] 


ql 


11     a 
A-j  i-rf  2-f  I  -exp[-crA]]-.5-[t  -«tp[-er2-A]])fl  =  —'[  .5  -exp[  -a-  A]  -{  .5-«p[  or  2-  A]]i3   =  .5-— •[  1  -  exp[   a- 2-  Aj' 
a  "  _/  or 


vl„V 


n  -H 
d  -H 


|ivln42-»2-v2nH{^2]2-v3n^l]j[^3]2-v3n^3]  jvln+2-^2-v2n^^2]2-v3n^lj-r3-[[v2n^2-v3n]-^3^2:2i: 
[vln42-^2-v2n-lt*2]2-v3ll-f<,l  4Tl]-[:*3l2-v3n-K,3-fr3J-[[v2n-r*2-v3n]-^3-K,2]2 


v2n-r^2-v3n]-^3-KI2]Tl 


vln+2-^2-v2n-(f^2]2-v3n-h,l  -tTl]-[[^3]2-v3n-K,3-rr3]-[[v2n-^2-v3n]-^3-K12J2 


r3Ii[v2°^2'v3°]^3^2]Mvi°+2^2,v2°^^2]2'v3°^iH^3]2'v3°'Kp]'>Ti'[^3]2'v3°^3]l 


vln42-02-v2o-lf*2f-v3n^l  -frIj-[[*3]   -v3n  -h,3  -rr3  ]  -[[  v2Q  -*2-  v3n]-*3  -K,2 


kln    =  V"          *    =  '  -5           W* 

=      \        -,0° 

Sensor  AOU  =  1  nm 

A  =  1-10~3       ■           A  =0.01      o 

A  =0.05     ■ 

A  =0.1 

0 

A  =0.5     ■ 

A  =1 

DifT|             ■                   Diflj             ■ 

Difl,              . 

Difq 

■ 

Din; 

o          Diflq 

-99.980002 

-98.037619 

-66.575636 

-49.955032" 

-45.804266" 

-12.535059" 

-33.264773 

-27.592024 

-6.501015 

-24.912739 

-18.198738 

-6.359547 

-19.892447 

-12  55896 

-6.09602 

-16.539079" 

8932986 

-5.155953 

-14.(38292 

6  525805 

-3.855277 

-12. 332923 

4  903215 

-2.5665*1 

-10.924554 

-3  8009(31 

-1.510314 

-9794141 

-3.04773 

-0.786681 

-8.865931 

2.529217 

0  356929 

-8  089429 

1.167554 

-0.137247 

-7.459679 

I  90955? 

-0.042686 

6.861717 

1.718907 

0.010612 

6.367239 

-1.57099 

-0.003816 

-S.93l5l7 

-  1 .449334 

0.004216 

-5.54706" 

1.343167 

-0.005 

-5.202706 

1.245689 

-0.004674 

4.893016 

-1.152855 

0.003587 

-4.612839 

1.062508 

-0.0O2352 

-4.358007 

0.973753 

-4.12511 

0.886509 

-J. 911334 

0801183 

-3.714333 

-0.718437 

-3532138 

-0.639022 

33.154822 


7.528865 


9.446393 
7.0516O5 


3.604176 


1.284654 
0.298693 


0.035621 

0.004605 


-0.008734 


1.894728 


4.579219 


0.002639 


0.008118 


5.581487-10 


0  469314 


0.821639 


-0.031427 


8.463799-10 


2.082199-10 
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Sensor  AOU  =  lOnm 


A  =0.01 


A  =0.05 


A  =0.1      o 


A  =0.5      . 


A  =1 


A  =5 


Difl; 

0 

Din; 

-  0Q  979987 

-  09  500457 

-49.9550M 

-48.894589 

-J3.2&7W 

31.73098 

-24.912764 

-22.951261 

19.892292 

-17.54656$ 

-16.53969 

-13  849295 

-14.137557 

-11.U552J 

-12.331612 

-9.078576 

-16.912565 

-7  456049 

-9.791141 

-6.140237 

-A.A<si754 

-5.(VT5oi7 

-8.083836 

-4.192852 

-7.422426 

-3.465175 

-6.852556 

-2.860974 

-6.355922 

-2.358669 

-5.918804 

-1  94688 

-5.530716 

-i.5938l4 

5.18351 

-1.36588 

-4.876759 

- 1 .067456" 

-4.58733 

-0.870461 

-4.329672 

0.708084 

-4.092597 

-6.574577 

-3.87511  i 

-6.465689 

-3.(574291 

-6.375535 

-3.488188 

-6.362477 

Difl; 


98  023693 


45.862954 


I7393T58 


18.165175 


-12.43854 


-8.68271 
6.118373 


4.328969 


3.065327 

2.16778 


■i.MAMI 


1.674191 
0.751363 


-0.522938 


0  362028 
-0.24926 


6.176673 


0.116232 

0.078746 


0.053094 


Difl; 


Difq 


Din; 


65  885691 


12.837194 


-5.38817* 


1.947916 
1.356479 


0.483994 


0.132892 
0.028186 


■6.005372 

0.001726 


32.043382 


-6.917853 


3.044708 
0.490318 


6.027738 


-1  873284 

-6.749395 

-0.013501 

-1.586052- 

io"< 

1.635119- 

10"< 

Sensor  AOU  =  lOOnm 


A  =0.1      a 


A  =0.5 


Din; 

■ 

Difl; 

-99  979843 

-99  484887 

-49  955004 

-48.896967 

-33.266816 

-3i  738627 

-24.912377 

-22.923018 

-19.89094 

17  45li78 

-16.535544 

-13.674001 

-I4.l3l8lj 

-10.895156 

-11.322629 

-8.769084 

10.909682 

-7io2l59 

-9  774059 

5.774855 

-8840153 

-4.76748 

-8.057607 

3  843648 

-7.391589 

-3.141649 

-6.817245 

2.569615 

-631636 

-2.102692 

-5.875281 

l  721164 

-5.483573 

-1.409027 

5.13312 

1.153674 

-48i75i6 

-  0.944669 

-4.531635 

0.773565 

4.271329 

-  0.6J347 

4  033201 

6518756 

-3.814445 

0.424819 

-3.61272 

-6.347893 

-3  42606 

0.284898 

A  =1 


Difl; 


-97.923273 


-45.855825 


-27  64985^ 


18.018772 


12.647194 


8.113983 


■5.470751 
3.685715 


2.480428 


T(56775g 
1.126598 


-0.75261 


0.505316 


0.339214 
■6.227685 


0.152814 


6.162559 
0.068829 


0.046191 
0.030999 


A  =5     . 

Dirt, 


A  =10     a 


Difl; 


A  =50      a 

Difl; 


-65.624462 


14.77194 


3.136706 
0.559593 


-0.091712 


-36.291752 

-4.450593 

-0.269804 

-0.611582 

-5.467805-10  ~A 

-6.493503 

-0.141654 

-2.531692- 

10'4 

-4.543601- 

lO"7 

-8.154349- 

io-lfl 
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APPENDIX  D:   KALMAN  GAIN  VERSUS  INTERARRTVAL  TIME 


This  appendix  shows  the  relationship  between  the  Kalman  gain  and  the  interarrival  time.   This 
information  can  be  used  to  select  a  range  of  time  values  to  tabulate  preprocessed  variance  data. 


n    =  0..25 


A    -  10 


u    -  50 


a    =  .23 


4$0 


rl    =  .25-u' 


r3 


2a 


vIq    -  rl  v2q    =  0         v3q    =  r3 


42  -  -•[  1  -exp[-o*A]l       ^3   -  expf-orAj 
a 


q3    =  .5-  —  •[l-«p[-0'2-A]] 


a 


1    rr 


A-   -•[[2-[l-exp[-o-A]]-.5-[l  -  exp[ -a-2-A]]ll 


q2   =  ^=-[.5-exp[-crA]-lf  .5-exp[-cr2-A]l] 


n-H 

D  -H  ! 
n-H 


j  vl„  +2-*2-v2„  j  *2]2-v3„  -KH  ]•[[  *3f-v3n  -Kp]  j  vln  +2-*2-v2n  j  *2]2-v3n  -rql  ]-r3  -[[  v2„  -Hfr2-v30]-»3  -h,2  j 
#1_  +2-^2-v2n  -If  *2]2-v3n  -»t,l  -frl  ]•[  j  *3  f-v30  -htf  -fr3  ]  -[[  v2_  -l*2- v3l-*3  -h,2  f 


LI     ■ 


v2n-^2-v3n]-03-ttj2|-rl 


r3 


/ln  -r2-^2-v2n  -If  ^2]2-v3n  -h,l  -frl  ]•[[  -tff  -y\  -h|3  «]-[[  v2B  4*2-v3nj-*3  -h,2]2 


r3il[v2°^2,v3°l^3^2]Mvl°+2^2'v2^^2l2'v3°^1H^3^v3n't<t3]-ffl,l^3]2,v3°'Kt3 


vl„  -*-*2-v2„  -If  42]  -v3„  -+ql  -h-1  ]•[[  *3  j   -v3n  -K,3  -rr3  J-[[  v2„  -*2-v3B]-eW  -h|2 


vl. 


kl 


v2n 

k2a       --  — ? 

■        r3 


n        rl 
KDATA    -  READPRN( Kalman) 

Sensor  AOU  =lnm 


v2n 

0        rl 


v3. 


k3n    =  ~y  APPENDPRN( Kalman)    =  [A    u     kl^     k2a25     k2b2J     k3-,5   . 


time    =  KDATA  <0>     Kl    -  KDATA  <> 


n    -  O..Ua{ 


KDATA 


o>l 
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KDATA        READPRN(Kalman5) 

Sensor  AOU  =5nm 


time    =  KDATA <0>     Kl     =  KDATA <2>  n    =  0  ..  last.  KDATA 


<0> 


K1n0.5  \- 


KDATA    =  READPRN(KalmanlO)  0>  <2  > 

'        lime    =  KDATA  ^  ^     Kl    =KDATA^ 

Sensor  AOU  =10nm 


n    =  Clast  KDATA 


<0>1 


K1n0 .5  U 


KDATA    --  READPRN(Kalman25)  0>  <s  > 

lime    =  KDATA  ^  *     Kl    -  KDATA  ^  * 

Sensor  AOU  =25nm 


n    =  0..u4  KDATA <0>] 


K1n0.5U 


tune. 
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KDATA        READPRN(Kalman50)  <0  >  <2  > 

time    -  KDATA  ^  *      Kl     =  KDATA    *  ' 

Sensor  AOU  =50nm 


n    =  0  ..la«t[  KDATA 


<0>" 


K1n0.5   - 


KDATA    =   READPRN(KalmlOO)  <>  <>  r  <Q>] 

v  '  time    =  KDATA ^^     Kl    =  KDATA ^  n    =  0  . .  1*4  KDATA ^ * \ 

Sensor  AOU  =  100nm 


K1n0.5  U 
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APPENDIX  E:  OFFSET  ERROR  OF  FILTER  POSITION  DISTRIBUTION 

This  appendix  provides  indications  of  how  far  from  the  target's  true  position  the 
filtered  position  is  displaced  by  the  IOU  process  adjustment  of  velocity  to  account  for 
the  possibility  of  a  manuever.   This  assumes  the  target  is  actually  moving  with 
constant  course  and  speed.   The  first  graph  cover's  the  entire  range  from  lnm  to 
lOOOnm  while  succeding  graphs  show  only  the  neighborhood  of  the  maximum  value 


Offset  Error  for  1 0  kt  Target. 

Interarrival  Time:  5  hours. 


10  100 

Sensor  AOU  2-sigma  radius  (in  nm). 


1000 


96 


Offset  Error  for  30  kt  Target. 

Interarrival  Time:  5  hours. 


.    40.4 


50 


55  60  65  70 

Sensor  AOU  2-sigma  radius  (in  nm). 


75 


Offset  Error  for  3  kt  Target. 

Interarrival  Time:  1  hour. 


Sensor  AOU  2-sigma  radius  (in  nm). 
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Offset  Error  for  1 0  kt  Target. 

Interarrival  Time:  1  hour. 


1.045 


11.5      12      12.5      13      13.5      14      14.5      15 
Sensor  AOU  2-sigma  radius  (in  nm). 


Offset  Error  for  30  kt  Target. 

Interarrival  Time:  1  hour. 


11       11.5      12      12.5      13      13.5      14 
Sensor  AOU  2-sigma  radius  (in  nm). 


14.5      15 
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Offset  Error  for  1 0  kt  Target. 

Interarrival  Time:  0.1  hour. 


0.0138 


0.0135 


1.05     1.1 


Sensor  AOU  2-sigma  radius  (in  nm). 


Offset  Error  for  30  kt  Target. 

Interarrival  Time:  0.1  hour. 


.  0.0414 

<  0.0413 

I  0.0412 

J  0.0411 

IB  0.041 
I 


0.0409 


Z  0.0408 
o 

w  0.0407 

8  0.0406 

<5 

°*  0.0405 


Sensor  AOU  2-sigma  radius  (in  nm). 
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APPENDIX  F:  ERROR  BETWEEN  FILTERED  AND  ESTIMATED  SPA  SIZE 

This  appendix  shows  the  difference  in  standard  deviation  between  the  actual 
filtered  position  distribution  and  the  estimated  distribution.  The  data  tabulated  shows 
more  complete  data  for  the  cases  modeled.   This  includes  the  observed  AOU 
(Obs.AOU)  and  the  predicted  AOU  (PRED  AO);  the  mean  value  of  the  standard 
deviation  calculated  for  30  runs  of  48  hours  of  300  observations  for  the  position 
predicted  by  the  IOU  model,  the  contacts  position,  and  the  filtered  position;  as  well  as 
the  actual  mean  speed  of  the  contact  over  the  30  runs  and  the  mean  filtered  value. 


4cn, 

Time= 

=4  Hours,  Target  Speed =10kts. 

Parameter*  Alpha-. 25.  V-IOMs. 

1  AT\. 

>' 

l«HJ 

^ 

^ 

lUU 
on. 

OU 

^^ 

40- 

~7^~ 

■*£s* 

20- 

I* 

0- 

i 

— i                 i                 i 

i 

100 


200  300  400 

Banaer  AOU  Standard  Damalion  w\  nm 


500 


600 


m. 


100 


30 


25 


I  20 


15 


10 


Time=1  Hour,  Target  Speed=3kts. 

Parwrtterv:  Aipna-25.  V-10  W» 


X* 

/' 

-^^ 

'                               jS^ 

X  "^^ 

^lx^^ 

/ 

1         i 1         i 1         i         i 

20 


40     60     80     100    120    140    160 

Saner  AOU  Standard  Davtttan  m  nm 


Time=1  Hour,  Target  Speed =10kts. 

Pararrwtors:  Alphas  25.  V-10  ktt. 


30 


25 


20 


15 


10 


y^^ 

X^ 

■***"               JF*^ 

x^ 

/ 

i                      i                      i                      I                      I                      i 1 

20  40 


60  80  100         120         140         160 

Sanaor  AOU  Standard  Daviabon  m  nm 
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Time=1  Hour,  Target  Speed =30kts. 

Alpha-  25.  V- 10  kt» 


40  60  80  100         120         140         160 

Sanaor  AOU  ttmttmi  Pwriton  w  m 


Time=.1  Hour,  Target  Speed =10kts. 

Paranmn:  Alpha- .25.  V- 10  kta. 
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Data  for  on  IOU  process  motion  model  beta  =25  and  sigma  based  on  10  knot  average  speed 
Time  interval  between  observations  of  4  hours: 


Standard  Deviation  of  error  in  nautical  miles: 

Speed  in  knots. 

Prediction  Contact 

Filter 

Contact 

Filtered 

mean 

mean 

mean 

SPA 

%Error 

Obs.  AOU 

Pred.  AO 

mean 

mean 

14.941 

11.957 

11.032 

11.499 

-4.06122 

25 

30.684 

2.959 

3.044 

20.041 

19.071 

16.716 

17.18 

-2.70081 

40 

34.216 

2.975 

3.145 

22.194 

23.148 

19.693 

20.594 

-4.37506 

50 

36.52 

3.078 

3.078 

33.145 

43.599 

32.113 

35.113 

-8.54384 

100 

47.234 

3.02 

2.635 

52.749 

81.301 

53.04 

59.109 

-10.2675 

200 

67.285 

2.92 

2.118 

92.721 

160.022 

93.447 

103.544 

-9.75141 

400 

108.788 

2.992 

1.881 

149.547 

245.976 

151.585 

148.352 

2.179276 

600 

152.841 

3.061 

1.895 

21 .938 

1 1 .925 

11.197 

11.499 

-2.62632 

25 

30.684 

10.133 

7.482 

27.577 

19.357 

17.561 

17.18 

2.217695 

40 

34.216 

10.026 

7.087 

27.27 

22.54 

19.703 

20.594 

-4.3265 

50 

36.52 

9.997 

6.81 

40.463 

43.917 

34.288 

35.113 

-2.34956 

100 

47.234 

9.952 

6.369 

57.402 

84.078 

54.921 

59.109 

-7.08522 

200 

67.285 

9.943 

5.999 

95.207 

154.341 

95.083 

103.544 

-8.17141 

400 

108.788 

10.307 

5.925 

146.205 

243.202 

147.258 

148.352 

-0.73744 

600 

15Z841 

10.076 

5.958 

56.41 

13.4 

13.682 

11.499 

18.98426 

25 

30.684 

29.933 

21.703 

59.391 

21.688 

22.186 

17.18 

29.13853 

40 

34.216 

29.931 

20.036 

63.26 

26.691 

26.537 

20.594 

28.85792 

50 

36.52 

29.962 

20.124 

78.533 

47.705 

47.533 

35.113 

35.37151 

100 

47.234 

30.019 

18.783 

101.683 

89.644 

78.28 

59.109 

32.4333 

200 

67.285 

30.001 

18.033 

133.696 

169.918 

121.254 

103.544 

17.10384 

400 

108.788 

29.908 

17.635 

169.943 

233.559 

160.386 

148.352 

8.111788 

600 

152.841 

30.072 

17.556 

Time  interval  between  observations  of  2  hours: 

Standard  Deviation  of 

error  in  nautical  miles: 

Speed  in  knots. 

Prediction  Contact 

Filter 

Contact 

Filtered 

mean 

mean 

mean 

SPA 

%Error 

Obs.  AOU 

Pred.  AO 

mean 

mean 

7.08 

4.859 

4.547 

4.658 

-2.383 

10 

14.219 

3.041 

3.374 

10.401 

9.202 

8.001 

8.471 

-5.54834 

20 

17.145 

3.017 

3.416 

16.528 

20.859 

15.156 

17.217 

-11.9707 

50 

24.218 

3.023 

2.904 

25.218 

40.585 

24.689 

28.317 

-12.8121 

100 

33.625 

3.068 

2.421 

43.418 

78.444 

43.725 

47.058 

-7.08275 

200 

50.636 

3.065 

2.144 

58.566 

115.538 

58.707 

64.879 

-9.51309 

300 

67.655 

2.897 

1.979 

9.704 

4.805 

4.583 

4.658 

-1.61013 

10 

14.219 

10.066 

8.383 

12.211 

8.83 

7.806 

8.471 

-7.85031 

20 

17.145 

9.972 

7.805 

19.499 

21.68 

16.535 

17.217 

-3.9612 

50 

24.218 

9.964 

7.004 

28.038 

39.535 

25.836 

28.317 

-8.76152 

100 

33.625 

10.025 

6.615 

43.761 

78.01 1 

42.761 

47.058 

-9.13128 

200 

50.636 

9.958 

6.384 

64.397 

115.971 

63.61 

64.879 

-1.95595 

300 

67.655 

9.975 

6.213 
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22.67 

5.196 

5.242 

4.658 

12.53757 

25.584 

10.066 

9.858 

8.471 

16.37351 

32.774 

22.442 

20.253 

17.217 

17.63373 

47.178 

43.851 

35.221 

28.317 

24.38111 

64.995 

83.895 

56.581 

47.058 

20.23673 

77.334 

119.128 

71.016 

64.879 

9.459147 

Time  interval  between  observations  of  1  hours: 
Standard  Deviation  of  error  in  nautical  miles: 


Prediction  Contact 

Filter 

mean 

mean 

mean 

SPA 

%Error 

3.607 

2.323 

2.146 

2.296 

-6.5331 

5.215 

4.445 

3.803 

4.159 

-8.55975 

8.437 

10.543 

7.439 

8.435 

-11.8079 

11.984 

19.878 

11.43 

13.716 

-16.6667 

18.691 

37.993 

18.477 

22.018 

-16.0823 

25.46 

58.34 

25.401 

29.352 

-13.4608 

4.568 

2.393 

2.243 

2.296 

-2.30836 

6.018 

4.461 

3.846 

4.159 

-7.52585 

9.107 

10.411 

7.472 

8.435 

-11.4167 

13.545 

19.994 

12  392 

13.716 

•9.65296 

19.93 

37.793 

19.213 

22.018 

-12.7396 

25.759 

54.988 

25.225 

29.352 

-14.0604 

9.34 

2.448 

2.452 

2.296 

6.794425 

10.781 

4.784 

4.551 

4.159 

9.425343 

14.764 

10.961 

9.315 

8.435 

10.43272 

20.459 

21.364 

15.801 

13.716 

15.20122 

28.395 

39.358 

24.728 

22.018 

12.30811 

36.583 

57.941 

33.572 

29.352 

14.37721 

Time  interval  between  observations  of  .5  hours: 
Standard  Deviation  of  error  in  nautical  miles: 


Prediction  Contact 

Filter 

mean 

mean 

mean 

SPA 

%Error 

1.03 

0.488 

0.474 

0.485 

-2.26804 

2.563 

2.203 

1.847 

2.033 

-9.14904 

3.622 

4.095 

3.014 

3.05 

-1.18033 

5.894 

9.608 

5.522 

6.688 

-17.4342 

8.469 

18.734 

8.256 

10.542 

-21.6847 

11.295 

27.866 

11.156 

13.745 

-18.8359 

1.533 

0.487 

0.475 

0.485 

-2.06186 

2.839 

2.205 

1.874 

2.033 

-7.82095 

4 

4.182 

3.131 

3.5 

-10.5429 

6.342 

9.788 

5.743 

6.688 

-14.1298 

10 

14.219 

30.061 

24.285 

20 

17.145 

29.933 

21.955 

50 

24.218 

30.047 

20.111 

100 

33.625 

29.957 

19.565 

200 

50.636 

29.938 

18.814 

300 

67.655 

29.979 

18.895 

• 

Speed  in  knots. 

Contact 

Filtered 

AOU 

Pred.  AO 

mean 

mean 

5 

6.429 

3.079 

3.489 

10 

8.15 

2.989 

3.515 

25 

12 

2.96 

3.119 

50 

16.684 

3.032 

2.737 

100 

24.227 

2.984 

2.329 

150 

31.113 

2.991 

2.215 

5 

6.429 

9.988 

8.922 

10 

8.15 

9.904 

8.387 

25 

12 

10.002 

7.7 

50 

16.684 

10.005 

7.255 

100 

24.227 

10.091 

7.009 

150 

31.113 

9.946 

6.891 

5 

6.429 

30.02 

26.103 

10 

8.15 

30.055 

24.381 

25 

12 

29.998 

22.467 

50 

16.684 

29.948 

21.356 

100 

24.227 

29.953 

20.68 

150 

31.113 

30.002 

20.393 

Speed  in  knots. 

Contact 

Filtered 

AOU 

Pred.  AO 

mean 

mean 

1 

2.162 

3.046 

3.178 

5 

3.755 

2.988 

3.437 

10 

5.499 

3 

3.298 

25 

8.187 

3.033 

2.875 

50 

11.759 

3.028 

2.58 

75 

14.768 

3.01 

2.435 

1 

2.162 

10.011 

9.607 

5 

3.755 

9.984 

8.905 

10 

5.199 

9.944 

8.461 

25 

8.187 

10.024 

7.945 

104 


9.101 

18.558 

8.7 

10.542 

-17.473 

1 1 .972 

27.983 

11.653 

13.745 

-15.2201 

3.513 

0.51 

0.513 

0.485 

5.773196 

4.761 

2.365 

2.207 

2.033 

8.55878 

6.158 

4.456 

3.791 

3.5 

8.314286 

9.329 

10.247 

7.273 

6.688 

8.74701 

13.198 

19.553 

11.421 

10.542 

8.338076 

16.49 

28.111 

14.945 

13.745 

8.730447 

Time  interval  between  observations  of  .1  hours: 
Standard  Deviation  of  error  in  nautical  miles: 


Prediction 

Contact 

Filter 

mean 

mean 

mean 

SPA 

%Error 

0.45 

0.43 

0.34 

0.377 

-9.81432 

1.155 

1.921 

1.059 

1.263 

-16.152 

1.689 

3.728 

1.621 

1.995 

-18.7469 

2.084 

5.502 

2.034 

2.573 

-20.9483 

0.495 

0.429 

0.348 

0.377 

-7.69231 

1.241 

1.934 

1.112 

1.263 

-11.9557 

1.761 

3.748 

1.665 

1.995 

-16.5414 

2.221 

5.549 

2.144 

2.573 

-16.6731 

0.778 

0.466 

0.428 

0.377 

13.52785 

1.721 

2.005 

1.368 

1.263 

8.313539 

2.408 

3.832 

2.104 

1.995 

5  463659 

2.947 

5.618 

2.685 

2.573 

4.352895 

Time  interval  between  observations  of  .05  hours: 
Standard  Deviation  of  error  in  nautical  miles: 
Prediction  Contact      Filter 
mean  mean  mean         SPA  %Error 

0.322  0.4  0.274  0.312     -12.1795 

83  1.843  0.795  0.978     -18.7117 


0.359 

0.41 

0.291 

0.312 

-6.73077 

0.871 

1.855 

0.82 

0.978 

-16.1554 

0.534 

0.442 

0.37 

0.312 

18.58974 

1.176 

1.908 

1.045 

0.978 

6.850716 

50 

11.759 

10.05 

7.625 

75 

14.768 

9.981 

7.512 

1 

2.162 

30.01 

28.592 

5 

3.755 

30.038 

26.093 

10 

5.199 

30.015 

24.905 

25 

8.187 

29.992 

23.28 

50 

11.759 

29.652 

22.652 

75 

14.768 

30.026 

22.335 

Speed  in  knots. 

Contact 

Filtered 

Obs.  AOU  Pred.  AO 

mean 

mean 

1 

0.592 

3.002 

3.192 

5 

1.504 

3.021 

3.081 

10 

2.22 

3.003 

2.847 

15 

Z78 

2.995 

2.734 

1 

0.592 

10.003 

9.568 

5 

1.504 

10.005 

8.966 

10 

2.22 

10.026 

8.711 

15 

2.78 

10.015 

8.563 

1 

0.592 

29.986 

28.414 

5 

1.504 

29.998 

26.532 

10 

2.22 

30.037 

25.88 

15 

2.78 

30.008 

25.47 

Speed  in  knots. 

Contact 

Filtered 

Obs.  AOU 

Pred.  AO 

mean 

mean 

1 

0.407 

2.971 

3.138 

5 

1.077 

2.975 

2.919 

1 

0.407 

9.996 

9.549 

5 

1.077 

10.028 

9.011 

1 

0.407 

29.997 

28.449 

5 

1.077 

30.011 

27.085 
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APPENDIX  G:   MODELS  OF  ACTUAL  ESTIMATED  FILTER  OPERATION 

This  file  simulates  the  basic  operation  of  the  ASSET  tracker.   A  contact  conducts  a  random  tour  at  constant  speed 
as  a  contact  would  in  the  simulation.   Time  interval  between  glimpses  is  an  exponentially  distributed  variable. 
Statistics  are  collected  to  see  how  well  the  tracker  tracks  the  target  over  the  course  of  its  path.   The  number  of 
points  plotted  was  adjusted  to  allow  details  to  be  seen. 

j    -  0..75         runtime    =  48  glimpse    -  .1  A.    =     ln(md(  1  ))•  glimpse    mi    -  4  '  v    -   10     uiiq    -  0   u    -   10  mg    -  u 

tm   M    -  im    -hAj  N     -  until  runtime  -  tin-  ,j  N2    -  last(  N  )  -  1   t    -  0  . .  N2        n    -  0  . .  N2  44 

Al     --     ln(rnd(l))-miA2    =     ln(md(  1  ))-rm  +A1       A3    -  -ln(rad(  1  ))-mi  4A2    A4    =  - ln( md(  1  )  )•  mi  +A3       A5    =  - ln( rnd(  I ) )•  mi  +A4 
A6    =     ln(md(l  ))-mi  4A5A7    =     ln(md(  I  ))-roi  +A6A8    =  -ln(md(  1  ))-nu  +A7A9    '-     ln(md(  1  ))-mi  +A8cl    =  2-  T-md(l  )c2    =2-T-md(l) 
c3    ■  2-T-md(l  )c4    --  2-T-rnd(l)c5    =  2-ir-md(l  )c6    =  2-  T-md(l  )c7    =  2-T-rnd(l  )c8    =  2-*-md(i)    c9    :2-T-md(l)    clO    -2-T-md(l) 
tc„    =  if  tmn  <&1  ,cl  .if  tm    <&2  ,c2  ,i(  un    <&3  ,c3  .iff  tm„  <&4  ,c4  .iff  un„  <A5  ,c5  .if  Unn  <16  ,c6  .iff  tmn  <Al  ,c7  .if  tm    <&%  ,c8  ,if|"  tmn  <&9  , 
randln    =  md(l)         ™«l2n    =  md(l)      rand3n    =  md(  1  )   r»nd4n    =  md(l) 

r  ,  _  l"  "    _  '         _  1  i 

xrngn    :     ^?   •    i     2-lni  rtndln     -co8i  2*TT«nd2n    yrngn    -     ^  ;•  ,{ -2*lnf  tmnd3n    >tini  2-T*  r»nd4n  i      rc»en    -  modj  tcn  -+  mdl  -     --j,2"Tj 

rspdn    -  v  -r(md(6)-3)       rbrgn    :2-T-md(l)     Yspdn    =  v«nfr  to,,1,"]        XspdQ    -  vcosi  ;  tcn  1 "  YditQ    =  0  XdiSQ    =0 

2 
100       ,         a 


Ydlsn  -H  '■ 

Yspd^+Ydis,, 

Xd"n  « 

Xspd^An+Xd..,, 

Y^n     ' 

V  '  S1TV    tC 

n  j 

Xspdn 

vcoi  tc_ 
n  : 

*n 

ymgn  +Yd..n 

lngn 

xrngn  +Xdisn 

y"Pdn 

rspdn-sini'rcseI1^| 

xrspd,, 

L"P<*a*C0*7pc,en]jj 

.    a    =    25   a   ^    ~       r3    =  —         rl    =  .25-uZ     vU    =  rl  v20    =  0   v30    =  r3     <f>2_    =  -•'■  1  -Mpf -flf  A_]] 
n  2  2*a  Of  l  *■  ■ ' 

*3n    ^expi-o-A,,:  qln   =  ^-j  A„-j  ±'[  [2«[  1  -«q»|  -d'A,]]-  .S«[  1  -exp[-a«2-  A,]]]] j  j 


2  2 

q2n   -  — •>  .5-expi-a-An;  -ff  5-expi j-CPl-A^]]]  q3Q    =  .5-  —  •!  1  -  expi    a-2-  A„ | 
a  a 


vln  +2-*2n-v2n  4  ^i  -v3n  -K,ln  ;•; ;  *3n     -v3n  +<,3n  +rJ  :-. ;  v2„  +*2n-v3fl  -*3n  +q2„l 

rl  * ' T^-~T '  Vr.      '   .2 '-^~. ~2 

;  vln  +2-^2n-v2n  4j  ^  .   -v3n  +qln  -rrl  -  •  *3„     -v3„  +q3n  4r3    -      v2„  ^2n-v3(1  •  *3n  +«,2nj 

v2n^2n-^n  -^n^ni-^-r 


„  _u       -  n^nn^n^n  ,-j  r  i 

n  ^  '  I  vl      -W-M?    -V?     4k  W.9         .„i      -Ml      4,1    !..      *%%      * 


v3n4H 


!  vln  +2-*2n-v2n  +  #2.]   -v3n  +qln  -frl  ■•   >3„     -v3n  ^3„  +r3  ■- ;  f  v2„  ^2n'v3n  ;.<*3n  +q2n  : 
!  vln  +2-^2n-v2n  +'  *2„  !2-v3n  -h,ln  -rrl  ,-!  ;  *3n  2-v3n  4q3      -      v20  ^2n-v3n  ^3„  +q2n;2 

r3'.= r^~r r^ — T ^ ' '~2 

vl0+2^2n-v2n  +  ^2„J   -v3n  ■*,!,, -frl    '  :^3n  ;   -v3n  4q3n  4r3    -      v2Q  4tf,2„-  v3n  ■  <^3„  4q2n 


«.>.4^,{°   ».'S^  Vyo-^    x0=ln8o  r0.^,    s;^{2) 
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i Y. 

41 

i x, 

+1 

lv*. 

41 

v*i 

-H 

Y  -a  4Vy,-  Vyexpi     CfA.    ]                                    "  Y  -a  4Vy   -  Vy  •  exp!  -a- A.  ,    '  . 
L-^~ =LLJ^,4r[««t41-Lj ~~- ^    4tt.t41-yr,pWVyt-«P[-aAJj 

X  ,a  4Vx,     ViL-exp     orA.1                                      X-a4V)L     Vt -expi  -a-A." 
— — £ 1 U  +kl,  4,  -jtag,  4,  -v ■ Li    +k2^41.    xnpdt41Vxt-exp.-a-Atj, 

|              r  Y,-a4Vyt-  Vy^exp'-a-A,":  • 
Vy.-exp.  -a- A,  ;  4k2b,  4,-  Ut,  4,  -I - = 4lc3,  41"   yrspd,  41  -  Vy^exp*  -a- A, 

[X  t-a4Vx,-  Vyexpf-a-A,'1" 

VXj-expr    a- A,'  4k2b,  M  •    tog,  ^ ' «    4fcJ,  ^  •    xnpd,  ^     Vyexp     OfA,    ! 


j  v'°"'t  -h  :      •  *»«  *'  *V  v2t  t  *2,  f-  v3,  4,1,  |  .  Xi0U{  ^  1      |  Xt  4i-[  1  -  exp:    a-  A,]  1-  Vx, 

Viou2,  +,      -!         ;  v2.  4tp2.-v3,-<p3,  4<j2,  !   v-                               i    ,. 

,+1;                 L      l        ^       tj^l^t  ;Yiout4H!       !  y, +!-•.  1-expl    a-A..   -Vy, 

1  ^-       i                    *V   -v3,  4q3,  :  V»out4H 

li         *        l                 j  i                                        exp  -a- A,   -Vx, 

exPr-a-At;Vyl 

Parameters  for  this  run: 

D    =  meao[A]    N2  =49      a  =7.071068  a  =0.25     u  =5 


AOU    -  Vmean(vl  )      IOU    -  ^mean(Vioul )     AOU  =1.087788  IOU  =1.197159 

Detennine  the  difference  between  the  true  target  position,  and  the  filtered  position: 

Fx,    =  X,  -  Xdu,  Fy,    =  Y,  -  Ydii, 

the  true  target  position,  and  the  contact  position  and  the  true  target  position,  and  the  IOU-predicted  position: 

Cx,    =  Ing,  -  Xdigj     Cy,    =  1m,  -  Ydi«,  Ix,    =  Xioiij  -  Xdii,  ly,    =  Yiou,  -  Ydi«, 

The  root-mean-square  error  is  found  for  the  data: 

F2x,    ^Fx,,2    F2y,   ^Fy,]2    C2x,    ^Cx,]2      C2y,    =  [Cy,j2      Ex,   =  rix,j2      Dy,   -  [lytf 
SFx    =  ZF2x       SFy    =  lF2y        SCx    =  ZC2x  SCy    -  lF2y  Six   =  II2x  Sly    =  II2y 


Fe,,    ,    |i^)4,(Sfy){..5       C«   =    f lggU(*gU J      tar    ■    (J^LhJSII.J        Ferr  =0.952507  lerr  =1.0.829 

^Uwt(Fx)    la*(Fy)i  4[  lut( Cx )    ta«( Cy )  j  j[lMt(Ix)    la*(ty)j 

Cerr  =1  740496 

The  mean  filtered  and  contact  velocities  is  found: 

FY,    ^  JiVx,;2  4f  Vy,"2|     CV(    --  j  xnpd,  J2  4f  ynpd,  j2       me*n(FV)  =8.574437  me*n(CV)  =10.102023 

mFV    -  mean(FV) 
*dev(FV)  =0.773034  *dev(CV)  =1.641369      mCV    -  mean(CV) 

The  RMS  error  in  filtered  position  is  compared  to  the  RMS  error  in  the 
contact  position: 

TR    =  Cetr  -  Ferr  TR  =0.787989 

The  difference  from  the  computed  error  is  then  computed  for 
the  current  run: 

COMPerr   -  AOU  -  Ferr        COMPerr  =0.135281 
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Filtered  track  using  the  actual  algorithm,  with  mean  glimpse  interval  0.05  hour  and  Sensor  AOU  =  5nm: 


Al  =0.8%146 

A2  =1.251034 

A3  =3.394831 

A4  =4.572997 

A5  =6.685208 

A6  =16.686341 

A7  =31.313094 

A8  =31.611555 

A9  =34.12239 


N2  =49 


4 


-2  0  2 

Xdis,  .tag,  ,X, 


xestQ    =  0 


y««o 


xenn       .  vln-  ,(   2-ln(md(l  )))-coa  2-T-md(l  )    +Xdi»n  ye«n    =  .  vln- -,(-2- ln(rnd(  I  )))•  sia  2- T-md(  1  )    +Ydisn 

xVe{    =  VXj-exp     a-A,    -Hl3,  ^-' xrspd,  ^  -  Vxj-exp;  -O- A,         yVe,    =  Vy,-exp   "CrA,    -Hd,  +,  •    yrspd,  ^  -  Vyt'expi     a*  A,  -I 
Vdiffx,    =  xVetVx,  mean(Vdiffx)  =0.133526  Vdifly,    =  yVe,  -  Vy(  mean(  Vdifly  )  =0.030425 


Filtered  track  using  estimation  technique  (using  full  noise  values): 

5 


Eatx.    -  xeal.  -  Xdia. 

Esty     =  yest.  -  Ydis. 

,2 
Ex,    =  ;  Eatx,       Ey,    =     Esty, 

SEx    =  ZEx         SEy    =  lEy 


-„    last(Ex)     laot(Ey) 
Eerr  =1  028216 
Ferr  =0.952507 


JnKan(vI  )  =1.087788 


-4 


-2  0  2 

Xdia, ,  log, ,  xeat, 
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Filtered  track  using  the  actual  algorithm,  with  mean  glimpse  interval  0.1  hours  and  Sensor  AOU  =  lOnm: 


o  5 


Al  =2.089507 
A2  =10.470267 
A3  =14  93814 
A4  =17  572047 
A5  =20.489058 
A6  =22.232939 
A7  =25.150709 
A8  =25.233842 
A9   =45  284693 


N2  =74 


B*0 


ye«o    =  0 


xestn    =  .jVln-;(-2-ln(rnd(l)))-co^2-T-™J(l)l+Xdisn  ye«n    =  ^vln^(-2-ln(nid(l  )))-st^  2-T-md(l  )  ]  4Ydisn 

xVet    =  Vjtj-exp     a-A,  !  -Hl3,  ^-i  xrspd,  ^  -  V^-expl  -a*  A,  ;  i     yVet    =  Vyt'expi  -or  A, ' -Hc3,  .^ -r  ynp^ +,  -  Vyt-expi  -a*  A,  '  ] 
VdifFx,    -xVej-Vx,  mean(Vdiffx)  =0.056692  Vdiffy,    =  yVe,  -  Vy,  me«n(  Vdiffy)  =0.013961 


Filtered  track  using  estimation  technique  (using  full  noise  values): 

80 


Eau.    -  xe«L  -  Xdit. 
E«tyt   =  yen,  -  Ydis, 

Ex,    =  [  Eax,  f   Ey,    =  [  E*yt 
SEx    =  ZEx         SEy    =  lEy 


-    ff    (SEx)        (SEy) 
""    "  4  Ust(Ex)    U«t(Ey)  ; 
Eerc  =1.925021 
Ferr  =2.165015 


*;me»n(vl)  =2.072638 


0  5 

Xdia,  .lngj.xet^ 
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Filtered  track  using  the  actual  algorithm,  with  mean  glimpse  interval  0.5  hours  and  Sensor  AOU  =  lOnm: 


£1  =9.779778 

A2  =19  746842 

A3  =23.257026 

A4  =23  336796 

A5  =25.985588 

A6  =26.130136 

A7  =  2"V  402486 

A8  =28.369044 

A9  =40.739811 


N2  =91 


150 
Xd«,,tagt,Xt 


300 


xeiiQ    =  0 


ye«o 


=  0 


xe»tn    =  .  vln--,(-2-ln(md(l  )))-co»  2-T-md(l)j +Xdisn  ye«n    =    !vla>  v(-2-ln(md(  1  )))-smi i  2-T-rnd(t  )    +Ydisn 

xVe,    =  Vxj-exp     a*  A,    +k3,  ^  •    xrspdj  +,  -  Vx^expf  -a*  A,  \]     yVet    =  Vy^exp1  -or  A,    "*3t  ^  •'  yrapd,  ^  -  Vy(-exp     a*  A, 
Vdiffx,    =xVe,Vxt  mean(Vdiffx)  =-0.104991  Vdifly,    =  yVet  -  Vyt  raean(Vdifly )  =0.040473 


Filtered  track  using  estimation  technique  (using  full  noise  values): 


Eatx.    =  xeat(  -  Xdia. 
Estyt    =  yestj  -  Ydiij 

Ex,    =     E**,f   Eyt    -r^f 
SEx    -  ZEx         SEy    =  lEy 


ton   -    Ti^),J^)1..5 
,    last(Ex)    la«(Ey)J 

Eerr  =3.402076 

Ferr  =2.944767 


;me*n(vl  )  =3.398984 


100  150 

Xdii|  .log,  .  xeat, 
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Filtered  track  using  the  actual  algorithm,  with  mean  glimpse  interval  0.5  hours  and  Sensor  AOU  =  30nm: 


Al  =3  915512 

A2  =691616 

A3  =15.287011 

A4  =20.200269 

A5  =30.204378 

A6  =34  696944 

A7  =35  833538 

A8  =38.376747 

A9  =38.669162 


N2  =96 


150 


Xdis, ,  lng,  ,  X( 


xe«0    =  0 


ye»0    =  0 


xe*n    =  .  vln-%(-2-ln(md(l)))-co»  2-T-md(l)i +Xdisn  yestn    =     vln- -/(-2-ln(rnd(  1  )))-sini  2-T-rnd(l  )    +Ydisn 

xVet    =  Vx,*exp     OfA,    -HL3,  ^  •    xrapd,  ^  -  VXj-expi  -a*  A,  ]  j     yVe,    =  Vyt-expi  -a*  A, 1  4k3,  ^  •    yripd,  ^  -  Vy,*exp     a-  A, 
VdiffXj    =xVe,-Vx,  mean(Vdiffx)  =-0.056681  Vdifty,    =  y Ve,  -  Vy,  mean(Vdiffy  )  =-0.145247 


Filtered  track  using  estimation  technique  (using  full  noise  values): 

50 


Eatx.    -  xest.  -  Xdis, 
Esty,    =  yest  -  Ydis, 

Ext    =     EstXj '      Ey,    =  ;  Esty, 
SEx    =  ZEx         SEy    =  ZEy 


.    lasi(Ex)     last(Ey) 
EerT  =7.774309 
Ferr  =5.817649 


/mean(vl  )  =7.596565 


Xdii,  ,  lag,  ,  xea. 
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Filtered  track  using  the  actual  algorithm,  with  mean  glimpse  interval  1  hour  and  Sensor  AOU  =  lOnm: 


300 


A I  =5.022727 

A2  =i  2.437775 

A3  =17  876472 

A4  =28  662312 

A5  =42.644194 

A6  =47.842995 

A7  =49.239743 

A8  =50.274019 

A9  =53  491883 


N2  =45 


-150  100 

Xdia,  .ln^  ,X, 


e«0    =  0 


ye«0    =  ° 


xestn    =  ,  vln-;(-2-ln(rnd(l)))-c^2-T-md(l)j +Xditn  yestn    =     vln'.,/(-2-  ln(md(l  )))•««  2-x-md(  1  )]  +Ydi»n 

xVe,    =  V^-exp   -o-A,    -*3t  ^-1  xnpdj  ^  -  Vjtj-expi  -a*A,    !     yVet    =  Vyt'expi  -a-  A,  I  -Hl3,  ^  ■!  yrspd,  +,  -  Vy^expi  -a*  A, 
Vdiffx,    =xVe,-Vxt  me«n(  Vdiffx )  =0.472115  Vdiffy,    =  yVe,  -  Vy,  mean(Vdifly  )  =-0.66914 


Filtered  track  using  estimation  technique  (using  full  noise  values): 


300 


Emx.    =  xe«t.  -  Xdi*. 
E*y(    =  ye«t  -  Ydia, 

Ex,    =     E*x,  2   Ey,    =  [  Eity, 
SEx    =  ZEx  SEy    =  ZEy 


c^   -    .T    (SEx)     ,  (SEy)  ".  " 
"       "    ,    last(Ex)    taM(Ey)! 

Eerr  =4.03991 

Ferr  =3.395647 


vmean(vl  )  =3.978631 


-150  100 

Xdit.  .Ins,  .  xest, 
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Filtered  track  using  the  actual  algorithm,  with  mean  glimpse  interval  1  hour  and  Sensor  AOU  =  30nm: 


Al  =0  155067 

A2  =4.15131 

A3  =9.030112 

A4  =14.972823 

A5  =19212 

A6  =21  285281 

A7  =24.284637 

A  8  =33.048225 

A9  =36.179972 


N2   =44 


50 


50  100 

Xdis,  .log.  ,X, 


200 


y««o 


--  0 


xestn    =     vln-  ,(-2-ln(md(l  )))-cosi  2-T'md(l  )    +Xdis 


yestn    =  ,  vln-v(-2-ln(md(t  )))-»ini  2-ttik1(1  )    4Ydisn 
xVe,    -  VXj-exp     a- A,    -Hc3,  +,  •    xrspd,  ^  -  Vx^exp     Of  A,         yVet    =  Vyt*exp!  -a*  A,  '  -HcJ,  +,  •    yrapd,  +,  -  Vy(-exp     orA, 
Vdiflfx,    =xVet-Vx,  mean(Vdiffx)  =-0.495172  Vdiffy,    =  yVe,  -  Vy,  mean(Vdifly )  =0.01928 


Filtered  track  using  estimation  technique  (using  full  noise  values): 

200 


Estx.    =  xe«,  -  Xdis, 
Estyt    =  yest,  -  Ydis, 

Ex,    =     Eax,       Ey,    =     Eaty, 
SEx    -  ZEx         SEy    -  ZEy 


Eert 


(SEx)    +(SEy)     .  5 


,    last  (Ex)    last(Ey)j 
Eerr  =9.81617 
Ferr  =12.186325 


mean(vl  )  =9.549774 


-50 


50  100 

Xdis,  .log,  ,xest( 


150 


200 
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Filtered  track  using  the  actual  algorithm,  with  mean  glimpse  interval  1  hour  and  Sensor  AOU  =  60nm: 


Al  =0.003175 
A2  =3.366652 
A3  =  11.418094 
A4  =14.422879 
A5  =19.245857 
A6  =24.664634 
A7  =26.370376 
A8  =30.333605 
A9  =30.916577 


N2  =54 


50  100 

Xdia,  .lag,  .X, 


xestr,     =   0 


ye*0 


xe»tn    =  ,  vl^-  ,(   2-ln(md(l)))-coa"2-T-nid(l):  +Xdisn  yesln    =  ,^- <<(-2- ln(ntd(  1  )))•  sin!  2-irmd(l  )'  +Ydisn 

xVe,    =  Vx,-exp     a*  A,     -HLJ,  ^  •    xrspd,  ^  -  Vx,-exp'  -a-  A,         yVet    -  Vy,-exp'  -a*  A,    -Hc3,  ^  •'  yrspd,  ^  -  Vy,-exp     or  A, 
Vdiffx,    =  xVe,     Vx,  mean(  Vdiffx)  =-0.202094  Vdiffy,    =  yVe,  -  Vy,  mean(Vdiffcr )  =0.294623 


Filtered  track  using  estimation  technique  (using  full  noise  values): 


Eatx,    =  xest.  -  Xdia, 
Estyt    =  yest,  -  Ydis, 

Ex,    =     Eatx,       Ey,    =  '  Eaty, 
SEx    =  ZEx         SEy    =  ZEy 


_    Een- 


(SEx)    ^(SEy) 


^last(Ex)    laat(Ey) 
Eerr  =15.785448 
Ferr  =12.430599 


mean(vl  )  =15.165514 


50  100 

Xdia, ,  tag, ,  xeat, 
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APPENDIX  H:  KALMAN  GAIN  VERSUS  SENSOR  AOU  SIZE 

This  appendix  shows  the  relationship  between  the  steady  state  Kalman  gain  Ka  and  the  Sensor  AOU  size. 
The  impact  of  using  the  the  approximations  from  Appendix  A  is  also  illustrated.   The  data  for  each  case  was 
written  to  a  seperate  file  and  the  graph  alone  displayed 
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Using  the  transition  and  process  noise  approximations  from  Appendix  A  to  simplify  the  calculations: 
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Kalman  gain  versus  sensor  AOU  in  nm  for  an  interarrival  time  of  0.01  hours: 
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Kalman  gain  versus  sensor  AOU  in  nm  for  an  interarrival  time  of  0. 1  hours: 
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Kalman  gain  versus  sensor  AOU  in  nm  for  an  interarrival  time  of  1  hour: 
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Kalman  gain  versus  sensor  AOU  in  nm  for  an  interarrival  time  of  2  hours: 
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Kalman  gain  versus  sensor  AOU  in  nm  for  an  interarrival  time  of  4  hours: 
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